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Abstract 

The Polynomial Chaos Expansion (PCE) technique recovers a finite second order random variable 
exploiting suitable linear combinations of orthogonal polynomials which are functions of a given stochas¬ 
tic quantity £, hence acting as a kind of random basis. The PCE methodology has been developed as 
a mathematically rigorous Uncertainty Quantification (UQ) method which aims at providing reliable 
numerical estimates for some uncertain physical quantities defining the dynamic of certain engineering 
models and their related simulations. 

In the present paper we exploit the PCE approach to analyze some equity and interest rate models 
considering, without loss of generality, the one dimensional case. In particular we will take into account 
those models which are based on the Geometric Brownian Motion (gBm), e.g. the Vasicek model, the 
CIR model, etc. We also provide several numerical applications and results which are discussed for a 
set of volatility values. The latter allows us to test the PCE technique on a quite large set of different 
scenarios, hence providing a rather complete and detailed investigation on PCE-approximation’s features 
and properties, such as the convergence of statistics, distribution and quantiles. 

Moreover we give results concerning both an efficiency and an accuracy study of our approach by 
comparing our outputs with the ones obtained adopting the Monte Carlo approach in its standard form 
as well as in its enhanced version. 


Introduction 

In this paper we present a Polynomial Chaos Expansion (PCE) approach to the solution of the following 
Ito stochastic differential equation (SDE) 

dX t = r(t, X t )dt + a(t, X t )dW t (1) 

where Xq is the initial value of the unknown stochastic process X t , the term r(t,X t ), resp. aft, X t ), rep¬ 
resents the drift, resp. the volatility of the process, while Wt is a standard Brownian motion and it is 
considered up to a certain finite time T > 0 to which we will refer, taking into account the financial setting, 
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as the maturity time, of, e.g., some underlying. According to ed Section 4.5], if the coefficients r(t,x) and 
a(t,x) are smooth enough, then we have the existence and uniqueness of the solution to ([I]). 

From a numerical point of view, between the most popular, and simple, method to numerically approx¬ 
imate the solutions for a large class of SDEs, there are the so-called Monte Carlo (MC) methods, which 
rely on pseudo-random sampling of independent increments of the Brownian Motion in 0 , see, e.g, |17j 
Chapter 9], f2T] Chapter 3] and references therein. 

The PCE technique is based on a radically different approach, namely it expresses the solution by means 
of polynomial basis of the probability space where the solution to 0 , at certain time T, is defined, see, e.g., 
ED] Chapter 3] for further details. 

Strictly speaking PCE recovers a random variable in terms of a linear combination of functionals whose 
entries are known random variables called germs or basic variables. Several method are available to com¬ 
pute these coefficients. In particular we choose the Non-Intrusive Spectral Projection (NISP) method which 
employs a set of deterministic realizations, see, e.g., |2U1 Chapter 3] for further details. Moreover the theory 
developed by Doss in |M], allows us to apply NISP avoiding unfeasible numerical problems of NISP approach 
to get numerical solution to (JT|). 

Originally the Polynomial Chaos idea was introduced by Norbert Wiener in his 1938 paper, see m, 
where he applies his generalized harmonic analysis, see, e.g., ED ■ Then P.D. Spanos and R.G. Ghanem in 
EH Section 2.4, Subsection 3.3.6], combined the PCE method with a finite element one, in order to compute 
the solution of PDEs in presence of uncertain parameters. Afterwards latter approach has been applied in 
several frameworks, e.g. to the simulation of probabilistic chemical reactions, see [5] , to stochastic optimal 
trajectories generation, see [4] , to sensitivity analysis, see [5], in order to refine specific numerical methods, 
see [B], and also to a rather large set of problems arising in engineering and computational fluid dynamics 
(CFD), see, e.g., [213 Chapter 6], [TyS] , and references therein. 

Our analysis has been focused on three well known one dimensional SDEs whose solutions are related, 
respectively, to the Geometric Brownian Motion equity model, see, e.g., uni Chapter 11, section 3], to the 
Vasicek model, see ED, and to the Cox, Ingersoll , Ross (CIR) interest rate model, see |T2]. We also refer 
to fig] Chapter 3, Section 3.2.3], for more details and references about the aforementioned models. 

To get an exhaustive comparison with already obtained results, e.g. by using Monte Carlo approaches, 
the PCE-approximation for each of the cited financial models, has been implemented for a set of three 
different volatility values, namely a = {15%, 25%, 30%}. Moreover we have analyzed the convergence of 
both the mean and the variance, to analytical values, as well as the distribution and two quantiles of the 
PCE-approximation, providing a wide test beach for PCE machinery. 

In order to show the advantages of the PCE approach, the results of our PCE-implementation have 
been compared versus the ones obtained using both the standard Monte Carlo (MC) and the quasi- Monte 
Carlo (QMC) techniques. The former is based on pseudo-random sampling and its convergence properties 
basically rely on the Law of large numbers and the Central limit theorem. The latter uses low-discrepancy 
sequences to simulate the process and it is theoretically based on the Koksma-Hlawka-inequality, see [9], 
which provides error bounds for computations involving QMC methods, see, e.g., ED Chapter 5]. 

We would like to underline that the first two cases, namely the gBm-model and the Vasicek model, are 
meaningful since they have an analytical solution to the related SDE, therefore we provide a rigorous study 
of the convergence property of the PCE method applied to them. In the CIR case, since we do not have an 
analytical solution, the application of the PCE technique allows us also to derive some considerations on 
the convergence properties of the related approximating solutions. 

The paper is organized as follows: in Section [T] the PCE approach is described in a general setting, while 
Section [2] deals with the NISP approach. Section [3] points out how to apply PCE in order to decompose 
the solution of the considered SDE, by using usual numerical method for solving ([!]), and how to employ 
Doss Theory in the PCE-approximation setting. In Section EDM we consider numerical applications of 
PCE technique, from the point of view of the Doss Theory, to approximate, respectively, the solutions to 
the gBm equity model, to the Vasicek model and to the CIR model. 

In section we give an overview of the PCE results in term of accuracy and computational time cost 
for each one of the considered models. 
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1 Polynomial Chaos Expansion 

Let (f2, £, P) be a probability space, where fl is the set of elementary events, £ is a cx-algebra of subsets of 
Q and P is a probability measure on £. 

Let us consider the Hilbert space of scalar real-valued random variables L 2 (fl, £, P), whose generic element 
is a random variable X defined on (f2, £, P), and such that 

E[X 2 ] = [ (X(w)) 2 dP(w) < +oo . 

Jn 

Notice that, as for Lebesgue spaces, the elements X £ L 2 (fl, £,P) are equivalent classes of random 
variables. Note that L 2 ( S1,£,P) is a Hilbert space endowed with the following scalar product 

E[XY] = (X,Y ) P = [ X(u])Y(u:)dP(uj) , 

J o 

being 

= E[^ 2 ] = J (.YH) 2 dPM , 

the usual norm. To shorten the notation ||X||p := ||X||^ 2 (q s p ^, and the related convergence will be always 
referred as mean square convergence or strong convergence. 


Among elements in L 2 (H,£,P) there is the class of basic random variables , which is used to decompose, 
as entries of functionals, the quantity of interest Y such as a random variable of interest, the solution of 
the SDE at time T. We notice that not all the functions £ : H —»• D can be used to perform such a 
decomposition since they have to satisfy, see, e.g., ns Section 3], at least the following two properties 

• £ has finite raw moments of all orders 

• the distribution function F^(x) := P (£ < x) of the basic random variables is continuous, _// being its 
probability density function (pdf). 

Since the increments of the Brownian motion are independent and normally distributed, from now on let 
us consider as basic random variable £ ~ Af( 0,1/2). The latter choice is motivate in Remark 1 of Sec. 

The elements in L 2 (H,£,P) can be gathered in two groups: in the first one we have the basic random 
variables , let us indicate them by £, ruling the decomposition, for simplicity let us call it Set £, while the 
second set is composed by generic elements, let us say Y, which we want to decompose using elements of 
the first set, which is referred as Sety. 

Let us denote by cr(£) the er-algebra generated by the basic random variable /, hence cr(£) C £. If we 
want to polynomially decompose the random variable Y in terms of £, then Y has to be, at least, measurable 
with respect to the cr-algebra <t(£). Exploiting the Doob-Dynkin Lemma, see, e.g., [23] Lemma 1.13], we 
have that Y is cr(£)-measurable, by detecting a Borel measurable function g : R —► R, such that Y = g(£). In 
what follows, without loss of generality, we restrict ourselves to consider the decomposition in L 2 (fl, <r(£), P), 
moreover the basic random variable £ determines the class of orthogonal polynomials {'Li(£)}ieN it is always 
indicated as the generalized polynomial chaos (gPC) basis. We underline that their orthogonality properties 
is detected by means of the measure induced by £ in the image space (D, 13(D)), where B(D) denotes the 
Borel a- algebra of D, in particular, for each i,j £ N, we have 

(^i>^j)p= / (€( u ))*j (£M)dP(w) = [ '& i (x)'& j (x)fz(x)dx. (2) 

Jn J n 

Since £ ~ A/"(0,1/2), the related set {Ti(x)}j g N is represented by the family of Hermite polynomials 
defined on the whole real line, namely D = R, and 
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(3) 






Figure 1: The Hermite polynomials up to degree 5 


Figure [T] provides the graph of the first six orthonormal polynomials, achieved by scaling each in © b y 
its norm in L 2 (fl, cr(£), P), namely, Vi £ N, Hq is divided by ||\Eq|| p := y/2 i i\. 

Latter polynomials constitute a maximal system in L 2 (fl, er(£), P), therefore every finite second order 
random variable Y can be approximated as follows 


= (4) 

2=0 

for suitable coefficients <q which depend on the random variable Y. We refer to eq. Q as the truncated PCE, 
at degree p, of Y. Exploiting previous definitions, taking i £ {0,... ,p}, and considering the orthogonality 
property of the polynomials we have 


:( y >^)p= TrrT2<ff> $ <)p> 


ll*i||p ll*iL 

and, since Y = g(£), we also obtain 

(Y, = (9, = [ g(x)'i/ l (x)f i (x)dx 

J q Jr 


(5) 

( 6 ) 


moreover Y bd converges in mean square sense to Y, see, e.g., [T5] Section 3.1]. The convergence rate of 
the PCE-approximation 0 in L 2 (fl, <r(£), P) norm is strictly linked to the magnitude of the coefficients of 
the decomposition. Indeed by the Parseval’s identity, we have 


+OO 


2=0 


furthermore, using the orthogonality property of Hermite polynomials in L 2 (f2, <r(£), P), the norm 
given by 


y(p) 


2 

p 


E^u 


2=0 


2 

P j 


of Q is 


then, exploiting the fundamental properties of the orthogonal projection in Hilbert Space, see, e.g., |351 
Theorem 4.11], we can estimate the mean square error as 


+ OO 


y-yW =\\y\\1- yw = E 


i=N+l 


(7) 


4 



























thus the coefficients rules the convergence rate. 

We would like to underline that the PCE of Y^ approximates the E-statistics using the Cj coefficients 
appearing in eq. 0- e.g. the first two centered moments are determined by 


E 



= Cq , 


( 8 ) 


Var 


y(p) 


p 

i= 1 




(9) 


2 Non-Intrusive Spectral Projection 

Let us consider the case of a functional M. which governs the dynamic of a quantity we want to study. We 
assume that M. is characterized by a random output Y depending on the input £, which is itself a random 
variable. The goal is to describe the behavior of Y in dependence of the variations of £, and with respect to 
the action of the functional At. In what follows, without loss of generality, we consider the 1 — dimensional 
case, indeed the results we shall obtain can be easily generalized to greater dimension. In our setting At 
will be the solution functional to an SDE describing a specific interest rate model , see later Sec. EDM 
for further details, and we also allow the output Y to depend on a set of parameters, let us call such a set 
0 £ M n , with n > 2, characterizing the interest rate dynamics, namely 

Y = M(£, 0) • 


In order to shorten the notation the parameters 0 are omitted in the At definition. Our aim is to exploit 
the so-called Non-Intrusive Spectral Projection (NISP) method to obtain the truncated PCE of the output 
y, taking into account that the basic random variable £ inherits all the assumptions made in the previous 
section. Therefore we have the following spectral projection 


y(p ) = I>*i(0, 

i =0 

where the coefficients are defined as in 0. 


We recall that, see, e.g. [2S1 Section 3.3], the NISP approach computes the scalar product in (§ by 
Gaussian Quadrature formulae, in particular, V* £ {0, ...,p}, we have the following Gauss-Hermite type 
result 


Ci 



N 


3 = 1 


( 10 ) 


where {CjlfLu { w j}jLi are the quadrature nodes and the weights of the one-dimensional Gaussian inte¬ 
gration rule. Moreover (101 requires the evaluation of the process At at a well defined set of realizations 
for the input random variable £. By the very definition of NISP, Y hh achieves spectral convergence with 
respect to the Y = Al(£) for increasing degree p, see, e.g., |2Dj Appendix B]. 

We note that the number of quadrature points N are not linked a priori to the degree p of the truncation. 
Instead, the precision of the coefficients Cj is related with N, thus a reasonable and standard choice that has 
used in literature, and which we will use in what follows, is to take N = p, since the polynomials involved 
are at least of degree p. Moreover, in what follows, we refer to {£j}jLi as quadrature nodes as well as a set 
of independent realizations of the basic random variable. 


2.1 Flowchart of the NISP computation 

In order to better explain the NISP approach, let us describe it exploiting the following flowchart 
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We would like to underline that the set of orthogonal polynomials, which constitute a basis in L 2 { fl, cr(£), P), 
is computed off-line once the basic random variable £ has been given and independently on the particular 
model we want to decompose. Analogously, the same happens for the set of N = p realizations of £ and 
weights needed to apply the Gaussian quadrature method. 

A different situation happens about simulation of the process and related computations. Indeed, we have 
to deal with data of the type {A^(£j)}^ x , which are sensible to changes of the model parameters. This 
is why such a difference has been highlighted in our flowchart using yellow, resp. blue, rectangles, which 
represent the off-line computations, resp. the model correlated computations. 


Remark 1 . Let us focus our attention of the basic random variable. Other possibilities are available, e.g. 
we can take £ ~ W(0,1), £ ~ exp(l), etc., nevertheless the Gaussian choice appears as the most natural one, 
since the Wiener increments are also distributed as Gaussian variables. In particular in our computations 
we consider £ ~ A/”(0,1/2), where the specification of the variance value is due to software restrictions. In 
particular we have used the Scilab toolbox which provides a powerful environment to implement PCE decom¬ 
position, also giving to the user a very flexible and efficient tool for the concrete computations involved by 
the PCE method itself, particularly by exploiting The Scilab toolbox, the Hermite polynomials require such 
input random variable in order to satisfy &> see F7R for further details. 

Remark 2 . The sampling of size L for the PCE ofYW) is achieved by detecting a unique sample {£/}/)=! 
of size L, of the basic random variable £, then we use eq. ( |10| |, to each realization, namely 

p 

Y i p) = $>*<(&) ,* = !> •••>£> 

2 — 0 
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therefore the collection 


N ri L ■ 


is the required PCE sampling. 


Polynomial Chaos Expansion and Stochastic Differential Equa¬ 
tions 


Let us consider a stochastic process {A t }t>o which satisfies the following SDE 

dX t = r(X t )dt + <j(X t )dW t , 


( 11 ) 


where Wt := {H't}te[o T] is a R-valued Brownian motion on the filtered probability space (fi, E, Ej,P), 
{S t } *e[o,T] being the filtration generated by W t . 


Our aim is to write the PCE approximation for the process AY, solution to (111, at a given, positive and 


finite time T, namely the random variable representing the process at that time. Without loss of generality, 


in what follows we consider the one dimensional version of the eq. (Ill, hence Xt is a scalar, real-valued 


random variable. In particular we will focus our attention on three specific models which can be described 


by the eq. (Ill, namely the geometric Brownian motion (gBm), the Vasicek model and the CIR interest 


rate model. Since in each of the aforementioned models the random variable Xt has finite variance, there 
exists a suitable probability space (fi, Ey^P), where Xt is well defined and such that Xt G L 2 (fl, Ey, P). 
Therefore the PCE method applies provided the definition of the the functional At such that Xt = M (£), 
for a suitable choice of the basic variable £ as seen in Sec. |2.1| 


3.1 Usual Numerical Methods for SDE and PCE 

Numerical methods used to approximate solutions to eq. 0, usually exploit the independent increments of 
the Wiener process. For instance, by considering the Euler-Maruyama scheme, see, e.g., [T71 Section 10.2], 
for each k = {0,..., L — 1}, we have 


Xt k+1 = A, 


tk 


r(tk, X tk )Atk + cr(tk, A t Jv / ^hA/'(0,1) 


where { t k } k=0 , L G N is a set of increasing time steps 0 = t 0 < ti,.. ., Y_i < t k = T, and A 0 is the value of 
the process at starting time t = 0, while At/ C = t k +i — tk- To shorten the notation let us define X k := X tk , 
for each time step, then the solution of (111 at time T reads as follows 


L—l 


Xt = X L =Tq (l + r(tk, X tk )At k + cr(£ fc , X tk )y/At k V2fkJ 


( 12 ) 


fc =0 


or equivalently 

X r=M(t) €=(&,..., &-i), 

where aforementioned set of i.i.d. Gaussian random variables with zero mean, and with 

variance equal to 1/2, namely they are independent A/"(0,1/2). Therefore we are in position to apply the 
PCE decomposition, see, e.g., El Secton 3.2] for further details concerning multivariate decomposition. 
We would like to underline that, even if we consider the Euler-Maruyama scheme, the previous approach 
can be applied considering any other method which is based on independent increments of the Brownian 
motion, such as Milstein method, see, e.g., [T71 Section 10.2, Section 10.3]. 

Let us recall that to get reliable results L has to be great enough, usually L = 100, which implies a high 
computational cost since the complexity of the NISP method increases dramatically when more than 20 
inputs are involved. The latter problem is referred as the curse of dimensionality , see |20| . 


3.2 Solution functional of a SDE 


The theory developed in ]33] gives conditions to express the solution of an autonomous SDEs as a functional 
of the Brownian motion. In what follows let us consider an autonomous stochastic differential equation (111, 
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by assuming that the drift and the volatility are Lipschitz real valued functions, the solution process {X t } t >o 


of (11) can be expressed as 


X t = H(D t ,W t ), te[0,T], 


where Wt denotes the standard Wiener process at time t and H : 
that 

( dH(x,y) _ 


^= a{H{x,y )) 
H(x , 0) = x 


(13) 

is a suitable function such 

(14) 


for every fixed x £ 


Moreover {D t } f > 0 hr (131 is a continuous process, adapted to the filtration of the Brownian motion, and 
such that for nearly all w £ 

D(u) = exp {- / 0 W * M a '(H(D(w),Wt(oj)))} L. Wt(u))) - W t (oj)))cr'(H(D( U ), Wt(u)))) ,t> 0 

D(lj) = Xq , t = 0 

(15) 

thus {D t } t >o is determined path-wise by means of the aforementioned deterministic ODE. Moreover the 
dependence on time t of D(ui) is skipped to shorten the notation. For further details on properties of these 
functions, see, e.g., [ 33] . 

Latter approach is suitable for the PCE-approximation method, since it defines the solution of the SDE 


(111 as a function of the Wiener process, see (131. 

In particular we can developed the following general road map , essentially based on Doss analysis, to 
apply the PCE-approximation: 


• Given the SDE (111, the solution functional H is detected by solving (141. 

• By the very definition of the Weiner process 

Wt~AT(0,t) ,te [0,T], 

therefore Wt = V / 2t£ 1 where the basic random variable £ is 7V(0,1/2), see Sec. [ 2 ] Exploiting such 
feature and integrating the differential equation (15), where the Brownian Motion is expressed as 
above, the stochastic process {D f } t >o is path-wise determined in terms of £. 

• Eventually, by means of the solution process is computed as Xt = H(Dt, Wt) for T > 0. The 
aforementioned steps allows to express Xt in terms if the basic random variable £, providing a NISP 
friendly definition of Xt, see Sec. [2] 

It is worth to mention that previous steps are discussed in full detail in what follows, for each numerical 
application considered, namely for each functional solution H(D t ,Wt ) related to the specific dynamics of 
the equity, resp. of the interest rate models, we have taken into account. 


4 PCE approximation of the Geometric Brownian motion 


In what follows we analyze the geometric Brownian motion defined as the solution of the following stochastic 
differential equation 

dSt = rStdt + aStdWt , (16) 

where r,a £ R + and Wt is the usual Wiener process. In particular let us focus our attention on St, for T > 0. 


The PCE method is based on the definition of a suitable process A4. Let us consider (161 and apply the 
approach described in Sec. |3.2[ where the solution at time T is defined as 

S T = H(D t , W T ) ■ (17) 

Since cr{y) = ay, by means of (14), H : R x R —► R is given by 


H(x, y) = xe' 


(18) 












Furthermore the random variable Dt is determined by integrating (151 up to time T. Thus for almost all 
u! £ 12, we have 


D(u) = exp {- aW t {u)}{rH{D{u),Wt{u})) - w))) 

D{u) = So 


t > 0 


t = 0 


(19) 


where D(u>) denotes the derivatives with respect to time of D t (uj), for every fixed sample event ui £ 12. 
We would like to underlying that the dependence on time, t, in previous equation is omitted to shorten the 
notation. Furthermore 

W t ~AT(0,t) ,t £ [0,T], (20) 

which allows to set Wt = x/2t£, and by plugging it into (191, the solution at time Dt can be expressed as 


a functional depending on the input random variable £, namely Dt = A4i(£, 0). 

We note that if the gBm-model is used to describe the behaviour of the interest rate S t , e.g. considering 
interest rate for equity markets, the related output at fixed time T, i.e. St, depends not only on £ but on 
other parameters, gathered by 0 = (r, a) £ R 2 , which characterize its dynamics. 


The solution at time T of (191 is numerically computed by means of an adaptive Runge-Kutta method 


of order 4 (RK4), where the absolute tolerance, resp. relative tolerance, is set as le-7, resp. le-5; see, e.g., 
[35J Section II. 1, Scetion II.4], for further details. 

Thus, by means of ((171, we obtain 


A4(£,0) ■=S T = D T e aWT =MiK,0)e^. 


( 21 ) 


To shorten the notation, in what follows the dependence on 0 is omitted. 


Turning back to the analysis of St, we apply the NISP method, exploiting eq. (10), in order to obtain 
the truncated following PCE 

5 T )= y>*i(o, 


3=0 


in particular let us set {<S't,.j}?L 1 = {■Ad^)}'^, where N = p, as the evaluation of the process (211 at 


defined by the Gaussian quadrature formula as stated in (101. Latter realizations are determined 
in two steps: 

• each Gaussian quadrature nodes, belonging to can be interpreted as the image of a suitable 

ojj £ fi through £. Therefore by integrating ( |19| ) for each path in {wi,...,Wjv} C f! we achieve a 
set of realization of Dt, let us say {A4i(^j)}v =1 . It is worth to mention that we are integrating N 
independent ODEs. 


By means of (211 we get the required set by point-wise multiplication of {A4i(^)}yL 1 by 


{< 




\ N 

1.7 = 1 


4.1 Numerical application: overview of the computations 


In what follows we give a PCE-approximation of the solution to the SDE (161, for some particular values 
of the parameters involved, see Table [T] 


Parameters r 

So 

T 

Values 3% 

100 

1 


Table 1: Parameters of gBm 

In order to make the discussion as complete as possible, our method is tested on a set of volatility values 

<r= 115%, 25%, 30%} . 
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Since eq. (16 1 has an analytical solution, such a solution will constitute our benchmark, in particular, at 

S T = S 0 e( r -^ T+aWT , 


time T, we have 


( 22 ) 


whose mean and variance are respectively 

E[S t ] =S 0 e rT , 

Var [S T ] = S$e 2rT (> 2t - l) . 

Then let us compute the absolute errors of the mean and the variance, namely 


Ap) 

t MEAN 

Ap) 

tyAR 


E [S r ] - E 


S, 


(p) 


Var [St] — Var 


Ap) 


the absolute value of the relative error being also considered. Due to the logarithmic scale, we look at their 
absolute values, highlighting the related order, rather than its numerical value. Therefore 


RE 


( p) 

MEAN 


RE, 


( p) 

VAR 


Ap) 

t MEAN 

E [Sr] 
Ap) 

t V AR 


Var [St] 

We note that such data have been computed for an increasing set of degrees p = {1, 2,..., 15} . 

For the aforementioned set of degree, and for each value of the volatility a, we will compute two quantiles 
Q 1 , where 7 = 99% and 7 = 99.9%, of the PCE approximation Sj?\ also writing their analytical values. 
Recalling that the standard statistics for Q 7 is the 7 -th sample quantile, namely the (K + l)-th realization 
of the sampling of Sj?\ sort in ascending order, such that K < [ 7 M], where M is the size of the sampling, 
and [■] denotes the integer part of the real number within the brackets. In our analysis we employ a Latin 
Hypercube Sampling (LHS) of S^\ see, e.g., j24]|25] and [26] for further references, of size M = 5000. Thus 
wg hcivc Kqq y 0 = 4951, while -^99.9% = 4996, see, e.g., m for further details, and we compute the absolute 
error of Q 7 versus the analytical values Q 7 . 

In particular we are aiming at computing the absolute error of Q 7 with respect to the analytical values, 
which are 


Q 1 = exp 


tVTz~ 


where Z 1 represents the quantile of the normal random variable of zero mean and unitary variance. Thus 
the error are defined as 

e 7 — Q '7 ; 

As last comparison let us estimate Q 7 , where 7 = 99% and 7 = 99.9%, by means of the aforementioned 


sample quantile applied, where a standard Monte Carlo sampling of the analytical solution St to (22 1 . 


The accuracy of such computation is represented by the standard error of the estimated average, namely 


for a set {QAf c (l)}i— 1 of L = 200 independent estimates of the quantile, whose arithmetic average is 


Q)f c = (1 /L) ■ 1 Qj C (0> we detect its standard error 


SE C 


VL 


where the estimated variance is 


a 


2 


1 

L- 1 


E(Qf c (0-Qf c ) 2 

1=1 
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Remark 3 . We would like to underline that the choice of the basic estimator for quantiles, i.e. the sample 
quantile, is due to focus our attention to the efficacy of the method, instead of taking care of the estimates 
accuracy as well as providing a fair comparison between the data achieved. Nevertheless the dedicated 
literature provides other techniques, which are often more accurate, such as the Two-phase quantile estimator 
presented in J38I/ . the L-estimator, or the d Harrel Davis (HD) estimators, see, e.g., for further details. 

4.2 a - 15% 

In this section we set er = 15%, while other parameters are displayed in Table [l] 

First let us display in Figure [2j Figure [3] and Table [2j the absolute and relative error of the average, resp. 
variance, of the PCE-approximation of the gBm. 




Figure 2: Semilogy scale plot of the absolute error of the mean (left) and the variance (right) computed via 
PCE for gBm at time T = 1, whose parameters are r = 3%, cr = 15% and starting value Sq = 100, for a set 
of degrees p = { 1,2,..., 15}. 
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Figure 3: Semilogy scale plot of the absolute value of the relative error of the mean (left) and the variance 
(right) computed via PCE for gBm at time T = 1, whose parameters are r = 3%, cr = 15% and starting 
value So = 100, for a set of degrees p = {1,2,..., 15}. 


Degree of PCE Average Error Variance Error Average relative error Variance relative error 


1 

4.3211e-03 

6.2663e+00 

4.1934e-05 

2.5934e-02 

2 

9.6288e-06 

6.2400e-02 

9.3443e-08 

2.5825e-04 

3 

7.8038e-08 

4.1041e-04 

7.5732e-10 

1.6986e-06 

4 

9.3644e-08 

1.6273e-06 

9.0877e-10 

6.7351e-09 

5 

9.3664e-08 

4.3078e-07 

9.0896e-10 

1.7829e-09 

6 

9.3664e-08 

4.3922e-07 

9.0896e-10 

1.8178e-09 

7 

9.3664e-08 

4.3925e-07 

9.0896e-10 

1.8179e-09 

8 

9.3664e-08 

4.3925e-07 

9.0896e-10 

1.8179e-09 

9 

9.3664e-08 

4.3925e-07 

9.0896e-10 

1.8179e-09 

10 

9.3664e-08 

4.3925e-07 

9.0896e-10 

1.8179e-09 

11 

9.3664e-08 

4.3925e-07 

9.0896e-10 

1.8179e-09 

12 

9.3663e-08 

4.3925e-07 

9.0895e-10 

1.8179e-09 

13 

9.3664e-08 

4.3925e-07 

9.0895e-10 

1.8179e-09 

14 

9.3664e-08 

4.3925e-07 

9.0895e-10 

1.8179e-09 

15 

9.3664e-08 

4.3925e-07 

9.0895e-10 

1.8179e-09 


Table 2: Absolute and relative errors of average and variance for PCE approximation of gBm at time T = 1, 
whose parameters are r = 3%, a = 15% and starting value Sq = 100 


Both average and variance error are stationary, despite an expected spectral convergence for increasing 
degree p. This can be motivated by the presence of an error, which is not due to the PCE method, that 
corrupts the expected spectral converge of . There are only two possible causes: the analytical approxi¬ 
mation of the solution St with ( |l7[ ) and the error coming form numerical computations of {A4(£.,)}^ =1 , see 
eq. (|2l}, used in (101. Let us focus our attention on the latter: as displayed in Figure [4] Figure [5] and Table 


I 


by increasing the precision used to compute {A / 1(^)}^ =1 the errors decrease. This is achieved by setting 
e absolute tolerance, resp. relative tolerance, of RK4 method, used to compute (191, at le-15, resp. le-10. 
Hence PCE errors are influenced by such approximation required to compute the coefficients. In Fig. [4] 
we can see the spectral convergence of the error for the PCE approximation up to 4-th degree for mean, resp. 
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up to the 6-th degree for variance. Note that, although the increased accuracy, the error displays again a 
stationary behavior, but actually no improvement on numerical methods can be implemented, which allows 
us to conclude that this is the only source of error. 




Figure 4: Semilogy scale plot for the absolute error of the mean (left) and variance (right) computed via 
PCE of gBm at time T = 1 ( whose parameters are r = 3%, a = 15% and starting value So = 100) with 
higher precision at computing {SrjjyLi 




Figure 5: Semilogy scale plot for the absolute value of the relative error of the mean (left) and variance 
(right) computed via PCE of gBm at time T = 1 ( whose parameters are r = 3%, cr = 15% and starting 
value So = 100) with higher precision at computing {Sr,j}jLi 
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Degree of PCE 

Average Error 

Variance Error 

Average relative error 

Variance relative error 

1 

4.3212e-03 

6.2663e+00 

4.1935e-05 

2.5934e-02 

2 

9.7228e-06 

6.2400e-02 

9.4355e-08 

2.5826e-04 

3 

1.5946e-08 

4.1085e-04 

1.5475e-10 

1.7004e-06 

4 

3.4008e-10 

2.0681e-06 

3.3003e-12 

8.5592e-09 

5 

3.2061e-10 

9.9706e-09 

3.1114e-12 

4.1265e-ll 

6 

3.2060e-10 

1.5312e-09 

3.1112e-12 

6.3373e-12 

7 

3.2060e-10 

1.5023e-09 

3.1112e-12 

6.2175e-12 

8 

3.2057e-10 

1.5025e-09 

3.1109e-12 

6.2186e-12 

9 

3.2057e-10 

1.5022e-09 

3.1109e-12 

6.2173e-12 

10 

3.2057e-10 

1.5020e-09 

3.1109e-12 

6.2165e-12 

11 

3.2048e-10 

1.5026e-09 

3.1101e-12 

6.2187e-12 

12 

3.2111e-10 

1.5039e-09 

3.1162e-12 

6.2241e-12 

13 

3.2078e-10 

1.5025e-09 

3.1130e-12 

6.2183e-12 

14 

3.2074e-10 

1.5029e-09 

3.1126e-12 

6 .2202e-12 

15 

3.2071e-10 

1.5019e-09 

3.1123e-12 

6.2160e-12 


Table 3: Absolute error of the average and the variance of PCE approximation of gBm at time T = 1 ( 
whose parameters are r = 3%, cr = 15% and starting value So = 100) for higher precision at computing 


Coming back to Table [lj let us compare a sampling of size 5000 of for p = 15, achieved with standard 
Monte Carlo technique, with the probability density function of the analytical solution to eq. (161 at time 
T, whose distribution is lognormal, namely 


X = Af ( log(5 0 ) + ( r ——IT 


ct 2 T 


St 


\ 


computed values are shown in Figure [6] were the standard Monte Carlo sampling of St is shown. 




Figure 6: Probability density function of gBm at time T = 1, whose parameters are r = 3%, a = 15% and 
starting value So = 100, (blue curve) and histogram of a Monte Carlo sampling (size = 5000) of the S^ 
for p = 15 (left). The right plot displays the analytical probability density function of gBm, for the same 
parameters values, and a standard Monte Carlo sampling of size 5000 of the gBm-analytical solution. 
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Let us determine the 99% and 99.9% quantiles of PCE-approximation, by means of the sample quantile 
based on Latin Hypercube Sampling (LHS) technique. These values are compared with the analytical ones, 
see Sec. 14.II for further details. 


Degree of PCE 

e 99% 

e 99.9% 

1 

5.6531e+00 

1.1050e+01 

2 

2.7440e-01 

6.0752e-01 

3 

8.1900e-02 

6.6493e-01 

4 

8.3246e-02 

7.5604e-01 

5 

8.1763e-02 

7.5891e-01 

6 

8.1672e-02 

7.5879e-01 

7 

8.1672e-02 

7.5878e-01 

8 

8.1673e-02 

7.5878e-01 

9 

8.1673e-02 

7.5878e-01 

10 

8.1673e-02 

7.5878e-01 

11 

8.1673e-02 

7.5878e-01 

12 

8.1673e-02 

7.5878e-01 

13 

8.1673e-02 

7.5878e-01 

14 

8.1673e-02 

7.5878e-01 

15 

8.1673e-02 

7.5878e-01 


Table 4: Absolute errors of the two quantiles Qgg% and <399.9% of the PCE approximation of gBm at time 
T = 1, whose parameters are r = 3%, a = 15% and starting value Sq = 100. 


Quantile absolute error 



Figure 7: Absolute errors of the two quantiles Qgg% and Qgg g % of the PCE approximation of gBm at time 
T = 1, whose parameters are r = 3%, a = 15% and starting value Sg = 100. 


Eventually let us compute the two aforementioned quantiles by means of a standard Monte Carlo sampling 
of the analytical solution St, see (22). As discussed in Sec. 
shown 


4.1 the standard errors of the quantiles are 


SE q mc = 0.0817117 

S'99% 

SE q mc = 0.2478969 

S'99.9% 
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4.3 <r = 25% 

In this section we set a = 25%, while the other parameters takes the values of Table [l] 


First let us display in Figure [8] Figure [9] and Table [5] the absolute and relative error of the average, resp. 
variance, of the PCE-approximation of gBm. 




Figure 8: Semilogy scale plot of the absolute error of the mean (left) and the variance (right) computed via 
PCE for gBm at time T = 1, whose parameters are r = 3%, a = 25% and starting value So = 100, for a set 
of degrees p = {1,2,..., 15}. 


Degree of PCE Average Error Variance Error Average relative error Variance relative error 


1 

3.2989e-02 

4.8289e+01 

3.2014e-04 

7.0513e-02 

2 

2.0607e-04 

1.3341e+00 

1.9998e-06 

1.9481e-03 

3 

8.2671e-07 

2.4397e-02 

8.0228e-09 

3.5625e-05 

4 

9.0548e-08 

3.3971e-04 

8.7872e-10 

4.9606e-07 

5 

9.3735e-08 

2.6355e-06 

9.0965e-10 

3.8484e-09 

6 

9.3744e-08 

1.2087e-06 

9.0973e-10 

1.7649e-09 

7 

9.3744e-08 

1.2457e-06 

9.0973e-10 

1.8190e-09 

8 

9.3744e-08 

1.2460e-06 

9.0973e-10 

1.8195e-09 

9 

9.3744e-08 

1.2460e-06 

9.0973e-10 

1.8195e-09 

10 

9.3744e-08 

1.2460e-06 

9.0973e-10 

1.8195e-09 

11 

9.3744e-08 

1.2460e-06 

9.0974e-10 

1.8195e-09 

12 

9.3744e-08 

1.2460e-06 

9.0973e-10 

1.8195e-09 

13 

9.3744e-08 

1.2460e-06 

9.0973e-10 

1.8195e-09 

14 

9.3744e-08 

1.2460e-06 

9.0973e-10 

1.8195e-09 

15 

9.3744e-08 

1.2460e-06 

9.0973e-10 

1.8195e-09 


Table 5: Absolute and relative errors of average and variance for PCE approximation of gBm at time T = 1, 
whose parameters are r = 3%, a = 25% and starting value S 0 = 100 


16 





































































Figure 9: Semilogy scale plot of the absolute value of the relative error of the mean (left) and the variance 
(right) computed via PCE for gBm at time T = 1, whose parameters are r = 3%, a = 25% and starting 
value So = 100, for a set of degrees p = {1,2,..., 15}. 


As in the case where cr = 15%, both average and variance error are stationary, despite an expected spec¬ 
tral convergence for increasing degree p. This can be motivated by the presence of an error, not due to PCE 
method, that corrupts the expected spectral converge of S^\ As before, there are only two possible causes: 
the approximation of the solution of fl!6|) w ith ( |17| and the error coming form numerical computations of 
{Ad^j)}.^, by means of (21 1 , used in (|lo|. Focusing the attention on the latter: as displayed in Table jbj 

by increasing the precision used to compute {Ad^)}^^ the errors decrease, namely we set the absolute and 
relative tolerance in RK4 method as le-15 and le-10. Hence PCE errors are influenced by such a approxima¬ 
tion required to compute the coefficients. In Tab. [6] we can see the spectral convergence of the error for the 
PCE approximation up to 4-th degree for mean, resp up to the 6-th degree for variance. Note that, although 
the increased accuracy, the error displays again a stationary behavior, but actually no improvement on nu¬ 
merical methods can be implemented, hence we can once again conclude that this is the only source of error. 
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Degree of PCE Average Error Variance Error Average relative error Variance relative error 


1 

3.2989e-02 

4.8289e+01 

3.2015e-04 

7.0513e-02 

2 

2.0617e-04 

1.3341e+00 

2.0007e-06 

1.9481e-03 

3 

9.2078e-07 

2.4398e-02 

8.9356e-09 

3.5627e-05 

4 

3.5168e-09 

3.4097e-04 

3.4129e-ll 

4.9789e-07 

5 

3.2968e-10 

3.8858e-06 

3.1993e-12 

5.6741e-09 

6 

3.2060e-10 

4.1608e-08 

3.1112e-12 

6.0758e-ll 

7 

3.2058e-10 

4.5733e-09 

3.1111e-12 

6.6780e-12 

8 

3.2057e-10 

4.2648e-09 

3.1109e-12 

6.2276e-12 

9 

3.2058e-10 

4.2643e-09 

3.1111e-12 

6.2268e-12 

10 

3.2055e-10 

4.2613e-09 

3.1108e-12 

6.2225e-12 

11 

3.2048e-10 

4.2659e-09 

3.1101e-12 

6.2291e-12 

12 

3.2108e-10 

4.2641e-09 

3.1159e-12 

6.2265e-12 

13 

3.2077e-10 

4.2605e-09 

3.1129e-12 

6.2213e-12 

14 

3.2074e-10 

4.2667e-09 

3.1126e-12 

6.2303e-12 

15 

3.2070e-10 

4.2598e-09 

3.1122e-12 

6.2203e-12 


Table 6: Absolute error of the average and the variance of PCE approximation of gBm at time T = 1, whose 
parameters are r = 3%, a = 25% and starting value So = 100, for higher precision at computing {Stj 


Coming back to Tablejlj let us compute a sampling of size 5000 of S'jP for p = 15, achieved with standard 
Monte Carlo technique, which is compared with the probability density function of the analytical solution 


to eq. (161 at time T, whose distribution is lognormal, namely 


X = AT ( log(S„) + ( r - ^ ) T , a 2 T ) , S T = 


x 


the computed values are shown in Figure [TO] were the standard Monte Carlo sampling of St is shown. 




Figure 10: Probability density function of gBm at time T = 1, whose parameters are r = 3%, cr = 25% and 
starting value So = 100, (blue curve) and histogram of a Monte Carlo sampling (size = 5000) of the S^’ 1 
for p = 15 (left). The right plot displays the analytical probability density function of gBm, for the same 
parameters values, and a standard Monte Carlo sampling of size 5000 of the gBm-analytical solution. 
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Then, for each degree p, we compute the two quantiles Q a for 7 = 99% and 7 = 99.9% by means of the 


sample quantile, see Sec. 4.1 based on Latin Hypercube Sampling (LHS) technique. Obtained values are 


shown in Figure ED and in Table 0 

Degree of PCE 


e 99% 


e 99.9% 


1 

1.6863e+01 

3.4495e+01 

2 

1.5613e+00 

5.1013e+00 

3 

1.7022e-01 

9.2564e-01 

4 

1.8926e-01 

1.6542e+00 

5 

1.7039e-01 

1.6940e+00 

6 

1.6839e-01 

1.6917e+00 

7 

1.6839e-01 

1.6911e+00 

8 

1.6840e-01 

1.6911e+00 

9 

1.6840e-01 

1.6911e+00 

10 

1.6840e-01 

1.6911e+00 

11 

1.6840e-01 

1.6911e+00 

12 

1.6840e-01 

1.6911e+00 

13 

1.6840e-01 

1.6911e+00 

14 

1.6840e-01 

1.6911e+00 

15 

1.6840e-01 

1.6911e+00 

quantiles Q 99 % and Q 99 . 9 % 

of the PCE approximation of gBm at time 


Quantile absolute error 



Figure 11: Absolute errors of the two quantiles Q 99 % and Q 99 . 9 % of the PCE approximation of gBm at time 
T = 1, whose parameters are r = 3%, a = 25% and starting value Sq = 100. 


Eventually let us compute the two aforementioned quantiles by means of a standard Monte Carlo sampling 
of the analytical solution of (16). The standard error of the quantile estimators are 


SE q mc = 0.1764685 

^99% 

SE 0 mc = 0.5685217 

^99.9% 
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4.4 <t = 30% 


In this section we set a = 30%, while the other parameters are taken form Table [I] 

First let us display in Figure [l2| Figure |T3| and Table [8] the absolute and relative error of the average, 
resp. variance, of the PCE-approximation of the gBm. 




Figure 12: Semilogy scale plot of the absolute error of the mean (left) and the variance (right) computed 
via PCE for gBm at time T = 1, whose parameters are r = 3%, a = 30% and starting value So = 100, for 
a set of degrees p = {1,2,, 15}. 


Degree of PCE Average Error Variance Error Average relative error Variance relative error 


1 

6.7908e-02 

1.0006e+02 

6.5901e-04 

1.0006e-01 

2 

6.1100e-04 

3.9769e+00 

5.9295e-06 

3.9770e-03 

3 

3.8351e-06 

1.0472e-01 

3.7217e-08 

1.0473e-04 

4 

7.4151e-08 

2.1059e-03 

7.1960e-10 

2.1059e-06 

5 

9.3718e-08 

3.2736e-05 

9.0948e-10 

3.2737e-08 

6 

9.3798e-08 

1.3416e-06 

9.1026e-10 

1.3417e-09 

7 

9.3798e-08 

1.8147e-06 

9.1026e-10 

1.8148e-09 

8 

9.3798e-08 

1.8204e-06 

9.1026e-10 

1.8205e-09 

9 

9.3798e-08 

1.8205e-06 

9.1026e-10 

1.8205e-09 

10 

9.3798e-08 

1.8205e-06 

9.1026e-10 

1.8205e-09 

11 

9.3798e-08 

1.8205e-06 

9.1026e-10 

1.8205e-09 

12 

9.3798e-08 

1.8205e-06 

9.1026e-10 

1.8205e-09 

13 

9.3798e-08 

1.8205e-06 

9.1026e-10 

1.8205e-09 

14 

9.3798e-08 

1.8205e-06 

9.1026e-10 

1.8205e-09 

15 

9.3798e-08 

1.8205e-06 

9.1026e-10 

1.8205e-09 


Table 8: Absolute and relative errors of average and variance for PCE approximation of gBm at time T = 1, 
whose parameters are r = 3%, a = 30% and starting value Sq = 100 
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Figure 13: Semilogy scale plot of the absolute value of the relative error of the mean (left) and the variance 
(right) computed via PCE for gBm at time T = 1, whose parameters are r = 3%, cr = 30% and starting 
value So = 100, for a set of degrees p = {1,2,..., 15}. 


As for the previous settings, i.e. a = 15% and a = 25%, both average and variance error are stationary. 
The motivations are the same: the analytical approximation of fl!6|) w ith (19) and the error coming form 
numerical computations of by means of (211, used inffiok. 

As shown in Table 9 the errors decrease by increasing the accuracy to compute the {A / 1(^ J )}^ =1 . Moreover 
we can see the spectral convergence of the error for the PCE approximation up to 4-th degree for mean, resp 
up to the 6-th degree for variance. Note that, although the increased accuracy, the error displays again a 
stationary behavior, but actually no improvement on numerical methods can be implemented, which allows 
to say that this is the only source of error. 
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Degree of PCE Average Error Variance Error Average relative error Variance relative error 


1 

6.7908e-02 

1.0006e+02 

6.5901e-04 

1.0006e-01 

2 

6.1110e-04 

3.9769e+00 

5.9304e-06 

3.9770e-03 

3 

3.9292e-06 

1.0473e-01 

3.8131e-08 

1.0473e-04 

4 

1.9968e-08 

2.1077e-03 

1.9378e-10 

2.1078e-06 

5 

4.0106e-10 

3.4562e-05 

3.8921e-12 

3.4563e-08 

6 

3.2092e-10 

4.8507e-07 

3.1144e-12 

4.8508e-10 

7 

3.2065e-10 

1.1972e-08 

3.1118e-12 

1.1972e-ll 

8 

3.2063e-10 

6.2851e-09 

3.1115e-12 

6.2852e-12 

9 

3.2064e-10 

6.2283e-09 

3.1116e-12 

6.2285e-12 

10 

3.2063e-10 

6.2219e-09 

3.1115e-12 

6 .2220e-12 

11 

3.2057e-10 

6.2321e-09 

3.1109e-12 

6.2322e-12 

12 

3.2115e-10 

6.2226e-09 

3.1166e-12 

6.2228e-12 

13 

3.2085e-10 

6.2201e-09 

3.1137e-12 

6.2203e-12 

14 

3.2078e-10 

6.2329e-09 

3.1130e-12 

6.2330e-12 

15 

3.2075e-10 

6.2165e-09 

3.1127e-12 

6.2167e-12 


Table 9: Absolute error of the average and the variance of PCE approximation of gBm at time T = 1 ( 
whose parameters are r = 3%, a = 30% and starting value So = 100) for higher precision at computing 
! s r..i\j , 


Coming back to Table [lj a sampling of size 5000 of S^ for p = 15, achieved with standard Monte Carlo 
sampling, is compared with the probability density function of the analytical solution to eq. (161 at time 
T, whose distribution is lognormal, namely 


X=Af log(S 0 )+ r-- T 


cr 2 r 


, St = e 


x 


the computed values are shown in Figure [14] were the standard Monte Carlo sampling of St is shown. 




Figure 14: Probability density function of gBm at time T = 1, whose parameters are r = 3%, a = 30% and 
starting value So = 100, (blue curve) and histogram of a Monte Carlo sampling (size = 5000) of the S'jP 
for p = 15 (left). The right plot displays the analytical probability density function of gBm, for the same 
parameters values, and a standard Monte Carlo sampling of size 5000 of the gBm-analytical solution. 
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Then, for each degree p, we compute the two quantiles Q a for 7 = 99% and 7 = 99.9% by means of the 
sample quantile, see Sec. 4.1 based on Latin Hypercube Sampling (LHS) technique. Obtained values are 
shown in Figure [15] and in Table [TO] 


Degree of PCE 

e 99% 

e 99.9% 

1 

2.5084e+01 

5.2332e+01 

2 

2.8318e+00 

9.8090e+00 

3 

2.2776e-01 

7.0824e-01 

4 

2.7654e-01 

2.2441e+00 

5 

2.2992e-01 

2.3467e+00 

6 

2.2387e-01 

2.3402e+00 

7 

2.2386e-01 

2.3380e+00 

8 

2.2393e-01 

2.3378e+00 

9 

2.2393e-01 

2.3378e+00 

10 

2.2393e-01 

2.3378e+00 

11 

2.2393e-01 

2.3378e+00 

12 

2.2393e-01 

2.3378e+00 

13 

2.2393e-01 

2.3378e+00 

14 

2.2393e-01 

2.3378e+00 

15 

2.2393e-01 

2.3378e+00 


Table 10: Absolute errors of the two quantiles Q gg% and < 599 . 9 % of the PCE approximation of gBm at time 
T = 1, whose parameters are r = 3%, a = 30% and starting value So = 100. 


Quantile absolute error 



Figure 15: Absolute errors of the two quantiles <5gg% and < 5 gg. 9 % of the PCE approximation of gBm at time 
T = 1, whose parameters are r = 3%, er = 30% and starting value So = 100. 


As described in Sec. 4.1 let us compute the standard error of the 99% and 99.9% quantiles, by means 


of L = 100 independent estimates. 

The achieved values are SE n Mc = 0.2235454, resp., SE n Mc = 0.7501973. 

^99% '<99.9% 
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5 Vasicek interest rate model 


In what follows we consider the Vasicek model, see m, and |161 Section 4.4.3] for further details, which is 
defined by 

dR t = {a- /3R t ) dt + adW t , (23) 

where a, [3 and a are positive, real-valued constants and Rq is the starting value of the unknown process. 


The analytical solution to the SDE (23), at time T > 0, reads as follows 


Rt = e /3T i?o + w (l — e /3T ) + 


0 -PT 


e ps dW(s) , 


P J o 

therefore Rt is normally distributed and its mean, resp. its variance, is given by 


(24) 


E [.Rt] = Roe P T + -(1 - e /3T ) , resp.by Var [. R T ] = — (1 - e 2,JT ) . 
P P 


Note that we have used the notation {R*}t>o ' n order to highlight that the Vasicek model (23) is used to 


describe the evolution of interest rates. We also underline that such a model belongs to the class of the 
so-called one factor short rate model which includes, e.g., the Rendleman-Bartter model, the CIR model, 
the Hull-White model, the Ho-Lee model, etc., since Rt depends only on one market risk factor, see, e.g., 
|28l 129] and da Section 3.2.1] . 

It is worth to mention that the Vasicek model solution allows negative values for the interest rate evaluated 
at given time T, hence its use in credit market framework is often criticized by practitioners. 


5.1 PCE approximation 

Analogously to what we have done concerning the study of the gBm model, see Sec. [4] we can rewrite the 
solution to (231 using a functional H such that 


Rt = H(D t , W t ) , (25) 

for every fixed time T, which is often referred as maturity time. In particular, exploiting the Doss approach 


described in Section 3.2 the scalar function function H(x,y) is determined by solving the ODE (221, hence 
obtaining 

H(x,y) = cry + x. 

Then the process {D t }t>o is path-wise determined by solving 


D{u) =a-PH(D(u)),W t (u)) t> 0 
D(ui) = Rq t = 0 


(26) 


for almost all w € D, where D represents the derivative with respect to time. Moreover such time dependence 
of D is omitted to shorten the notation. 

As in Sec. [4] since the Wiener process is 


W t ~M{0,t) ,t G [0,T], 


(27) 


then we set Wt = V / 2and we consider the basic random variable £ to be Gaussian, namely we take 
£ ~ A/”(0,1/2) in order to satisfy the restrictions imposed by the software we have used to numerically 
implement the PCE method. 


Proceeding as we made studying the gBm-model, the solution at time T of (261 is numerically computed 


by means of an adaptive Runge-Kutta method of order 4 (RK4), where the absolute tolerance is set as le-7 
and the relative one is le-5, see, e.g., Ea Section II.1, Scetion II.4], for further details. Even if an analytical 
solution is available of (261, we use the numerical one in order to be coherent with the discussion of Sec. |4] 
Then by plugging (27) into (26), and integrating it up to time T, the solution Dt can be expressed 


as a functional depending on the input random variable £, namely Dt = A4i(£, 0), the latter input 
0 = (a, /?, a) £ R 3 collects the parameters that characterizes the dynamics of the Vasicek interest rate 
model. Due to (271, the eq. p5| becomes 


M(f, 0) := R t = aW T + A4i(£, 0) = ay/2Tf + M^, 0) . 


(28) 
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As made before, in what follows we shorten the notation by omitting the dependence on vector parameters 0. 


Applying the NISP procedure up to degree p, the PCE decomposition of Rt is given by 

v 


R, 


(. p) 


= 5>^(0, 


i —0 


where the coefficients are detecting using (101 and evaluating Rt in (28) on a set of realizations of the basic 
random variable £, namely considering N values Moreover, in what follows, we set N = p. The 

obtained simulated values for Rt will be indicated as follows 


{Rt^U = {Mm? i 


N 


The aforementioned realizations of the (281 are determined in two steps: 

• each Gaussian quadrature nodes, belonging to can be interpreted as the image of a suitable 

ujj £ through £. Therefore by integrating (261 for each path in {wi,... ,w/v} C Q we achieve a 
set of realization of Dt , let us say {AllIt is worth to mentioned that we are integrating N 
independent ODEs. 

• By means of (28) we get the required set {Rr,j}f— i, where the first summand {T\/2T£j}:?Li- 


5.2 Numerical applications: overview of the computations 

Let us implement the theoretical approach developed in previous section by a working example which 
consists in considering the Vasicek model for a set of volatility 

o- = {15%, 25%, 30%} , 

while the other parameters as set according to Table ED 


Parameters 

a ft 

Ro 

T 

Values 

0.1 0.2 

110 

1 


Table 11: Parameters and initial value for the Vasicek interest rate model. 


First we compute the absolute error as well as the relative error of the mean and variance of R ,?' 1 , with 
respect to degree p = {0,1,2,..., 15}. In particular, for every given value p, the absolute errors for the 
mean and for the variance are as follows 


Ap) _ 

t MEAN ~ 


IE [Rt] — IE 


R 


(p) 

T 


Ap) _ 

t VAR — 


Var [Rt] — Var 



while the relative errors are given by 


RE 


(p) 

MEAN 


RE 


(p) 

VAR 


ftp) 

t MEAN 

E [-Rt] 


e VAR. 

Var [Rt] 


Due to the logarithmic scale employed, we consider their absolute values, highlighting the related order, 
rather than their numerical value. 


Then let us compute the two quantiles Q 7 , where 7 = 99% and 7 = 99.9%, of the PCE approximation 
Rt\ where p = 0,1,...,15. Proceeding as in Sec. [4J we consider the 7-th sample quantile, namely the 
(I\ + l)-th realization of the sampling of R ^ such that K < [7 M], where M is the size of the sampling, 
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and [■] denotes the integer part of the real number within the brackets. In what follows we always employ a 
Latin Hypercube Sampling (LHS) technique of R^\ see (24, :25| and also [26] for further references, of size 
M = 5000. In particular we have K 99 % = 4951, while K 99 9 % = 4996, and we compute the absolute error 


of Q 7 versus the analytical values Q 1 . Therefore we evaluate e 7 = 


Q'y Q~) 


where Q 7 is the quantile 


of a normal random variable of mean E [Rt] and variance Var [Rt] and, indicating its cumulative density 
function by Fr 7 , , we have Q 1 = F^ : (7). 


Furthermore the same quantiles are computed by means of standard Monte Carlo approach using the 


analytical solution Rt to (231. The accuracy of the achieved values is expressed in terms of their standard 


error SEqmc , for 7 = 99%, 99.9%. Therefore the sample quantiles are computed L = 200 times, getting 
{Qdf c (l) independent estimates. Then the standard error reads as 


SEr 


Vl 


where a 2 = ^ £f =1 (Of C ( 0 - Q MC ^ 


and Q 7 c is the arithmetic average of {Q 7 c (0}^=i- 


Remark 4 . We would like to note that the choice of the basic estimator for quantiles, i.e. the sample 
quantile, has been taken to focus our attention to the efficacy of the method instead of its estimates accuracy 
as well as to provide a fair comparison between the data achieved. Nevertheless the dedicated literature 
provides more accurate techniques, such as the Two-phase quantile estimator, presented in jSS'j . the L- 
estimator and Harrel Davis (HD) estimators, see, e.g., for further details. 


5.3 (T - 15 % 

In this section the value of the volatility is set to a = 15%, while the other parameters are chosen as in 
Table E] 


First let us display in Figure [16] the absolute and relative error of the variance, while in Table [12] are 
collected both the absolute and absolute value of relative error for mean and variance concerning the Vasicek 
interest rate model. 


Degree of PCE 

Average Error 

Variance Error 

Average relative error Variance relative error 

0 

8.2289e-08 

1.8544e-02 

9.1279e-10 

1.0000e+00 

1 

8.2303e-08 

1.2490e-03 

9.1295e-10 

6.7349e-02 

2 

8.2312e-08 

1.2490e-03 

9.1305e-10 

6.7349e-02 

3 

8.2348e-08 

1.2490e-03 

9.1345e-10 

6.7349e-02 

4 

8.2321e-08 

1.2490e-03 

9.1315e-10 

6.7349e-02 

5 

8.2349e-08 

1.2490e-03 

9.1346e-10 

6.7349e-02 

6 

8.2328e-08 

1.2490e-03 

9.1322e-10 

6.7349e-02 

7 

8.2308e-08 

1.2490e-03 

9.1300e-10 

6.7349e-02 

8 

8.2327e-08 

1.2490e-03 

9.1321e-10 

6.7349e-02 

9 

8.2300e-08 

1.2490e-03 

9.1291e-10 

6.7349e-02 

10 

8.2298e-08 

1.2490e-03 

9.1289e-10 

6.7349e-02 

11 

8.2352e-08 

1.2490e-03 

9.1349e-10 

6.7349e-02 

12 

8.2311e-08 

1.2490e-03 

9.1303e-10 

6.7349e-02 

13 

8.2364e-08 

1.2490e-03 

9.1362e-10 

6.7349e-02 

14 

8.2321e-08 

1.2490e-03 

9.1315e-10 

6.7349e-02 

15 

8.2357e-08 

1.2490e-03 

9.1355e-10 

6.7349e-02 

Table 12: Absolute error of the average and the variance of R^ at time T = 
model parameters a = 0.1, (3 = 0.2, a = 15%. 

1, with respect to the Vasicek 
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Figure 16: Semilogy scale plot of the absolute error of the variance of R^ evaluated at time T = 1 (left) 
and the absolute relative error of the variance (right) of at time T = 1. In both cases the parameters 
of the Vasicek interest rate model are a = 0.1, /3 = 0.2, cr = 0.15. 


The absolute and relative errors of mean and variance are stationary, see Table 12 This can be motivated 
by looking at the solution to (28). By integrating (26 1 , upon plugging the definition of H{x,y), we get 

RRt = ctW t + e" /3T i?(0) + ^ (1 - e -/3T ) - aPe~^ T f W s e ps ds , 

P Jo 

we call such solution as RRt , in order to distinguish it from the analytical solution of the Vasicek interest 
rate model (23). Then, exploiting ([27]) , W s = V2s£, and Wt = V^Tt; 


(29) 


RRt = ctWt + e /3T f?(0) + ^ (l — e ,3T ) — £ • ^er/?e ,3T J \/2se^ s ds^ . 


Notice that the same random variable £ is used both for Wt and W s definition, since the latter is involved 
by detection of Dt, which is computed separately from Wt, thus the former is not linked with the latter. 
Therefore we can gather £ and due to £ ~ A/”(0,1/2), we conclude that RRt is normally distributed with 


E [RR t ] = e~ pT R{ 0) + ^ (l - e~ 0T ) , 
Var [RRt\ = o 2 ^Vt — /?e _/3T J 


(30) 

(31) 


Let us consider the PCE-approximation R^ \ for a fixed degree p, then by triangular inequality 


E [-R 7 1 ] — ^ 


R, 


( p ) 


< 


E — E 


R 


(p) 


+ |E[f?i? T ]-E[f? T ]| , 


The last summand on the right is null, by means of (|30j), therefore the error of the mean is driven by the 
first summand. Due to (|8j) and (10), the numerical approximation of {Rtj} jLi influences the accuracy of 
PCE-statistics, therefore such error corrupts the spectral convergence of the PCE method. Indeed in Tab. 
[T2] it becomes stationary at a value close to absolute tolerance of the RK4 method, i.e. le-7, which is used 
to compute the {Rryl/Li • 
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The absolute error of the variance is 


,(p) 

t VAR 


Var [Rt] — Var 




(p) 


< | Var [Rt] — Var [-R-Rt]| 


Var 


R, 


(p) 


- Var [RR t ] 


(32) 


Due to spectral convergence of PCE-approximation, the second summand on the right hand side is irrelevant 
if compared with the second one. Therefore let us focus our attention on |Var [Rt] — Var in particular 

let us compute the Taylor expansion, with respect to fi and centered at fi = 0 , of Var [Rt], resp. of Var [RRt], 
so that the variance of (23) can be approximated by 


Var[flr] = g(l-e-^) 


2 /? 


1 - (1 - 2/3 T + 


4/3 2 T 2 


) + OQ 3 2 ) 


= a 2 T-a 2 T 2 /3 + 0(/3 2 ) , 


while the Taylor expansion of (31) reads as follow 


Var [RRt] = o 2 ^Vt — fie dT J y/se^ s ds S j 
= a 2 (VT - ?/3T 3 / 2 + 0(0 2 ij 


= a 2 [T--fiT 2 )+0(f] 2 ), 


therefore 


,(p) 

t VAR 


<j 2 T - <j 2 T 2 /3 - <t 2 T + -a 2 T 2 fi + 0{f3 2 ) 


< \t 2 p + o{p 2 ) 


Hence let us compute e~y\ R for a set of decreasing values of /3 

fi = 10 [_5:O - 5:_1] , 

while the other parameters are a = 15%, a = 0.1, Rq = 110 and T = 1. The results are displayed in Figure 


(33) 


17 and we note that they agree with the theoretical bound given by (33). 



Beta 


Figure 17: Log-Log scale plot of the variance absolute error of R^fjf 1 Fl> at T = 1 for different fi = 10^ 5:0 ' 5: ^ 
compared with order 1. The parameters of the Vasicek interest model are a = 0.1, er = 0.15 and R 0 = 110. 
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Remark 5 . By means of (17]) the error of the variance is meaningful for the convergence in mean square 

(p) '—' 

sense of Rf, to Rt as p —► +oo. Moreover such computations show that the PCE-approximation actually 


achieves spectral convergence, but to (29 1 , which differs in mean square norm to the analytical solution Rt, 
see eq. (241, of |Var [Rt] — Var [RRt] , hence motivating the behaviour of the error for variance in Fig. 16 


Remark 6 . We note that a does not influence the computations shown above, therefore the same procedure 
and related results are also valid for the other values of volatility, whose associated analyses are therefore 
skipped. 


Coming back to the values in Table 


11 


?(p) 


let us compute a Monte Carlo sampling of size 5000 of Rf, 
with p = 15, whose histogram, see fig. |18[ is compared with the probability density function of the normal 
random variable Rt- 

Moreover in the same fig. [18] we show a sampling obtained exploiting the standard Monte Carlo technique. 




Figure 18: Left plot: Probability density function of the Vasicek interest rate model, with parameters 
a = 0.1, /? = 0.2, a = 0.15 and Rq = 110, (blue curve) and histogram of a Monte Carlo sampling (size 
= 5000) of the Rj?) for p = 15. With the same parameters, in the right plot we compare the analytical 
probability density function of the Vasicek model, with the one obtained using the Monte Carlo sampling 
of size 5000, for Rt, being T = 1. 


The errors of quantiles for 7 = 99%, resp. for 7 = 99.9%, are shown in Fig. 19 and in Table 13 
values are estimated by means of the sample quantile, see Sec. |5.2|for further details. 


These 
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Degree of PCE 

e 99% 

e 99.9% 

0 

4.3141e-02 

5.7307e-02 

1 

2.6330e-01 

3.5319e-01 

2 

2.6330e-01 

3.5319e-01 

3 

2.6330e-01 

3.5319e-01 

4 

2.6330e-01 

3.5319e-01 

5 

2.6330e-01 

3.5319e-01 

6 

2.6330e-01 

3.5319e-01 

7 

2.6330e-01 

3.5319e-01 

8 

2.6330e-01 

3.5319e-01 

9 

2.6330e-01 

3.5319e-01 

10 

2.6330e-01 

3.5319e-01 

11 

2.6330e-01 

3.5319e-01 

12 

2.6330e-01 

3.5319e-01 

13 

2.6330e-01 

3.5319e-01 

14 

2.6330e-01 

3.5319e-01 

15 

2.6330e-01 

3.5319e-01 


Table 13: Absolute errors of the two quantiles Qgg% and Qgg.g% of the PCE approximation of the Vasicek 
Interest rate model at time T = 1. The parameters are set as a = 0.1, /3 = 0.2, a = 0.15 and Rg = 110. 


Quantile absolute error 



Figure 19: Absolute errors of the two quantiles Qgg% and Qgg.g% of the PCE approximation of the Vasicek 
interest rate model at time T = 1. The parameters are set as a = 0.1, (3 = 0.2, a = 0.15 and Rg = 110. 


The quantiles computed by means of a standard Monte Carlo sampling of the analytical solution, whose 
size is equal to the one used for PCE, give the following standard error 


SE q mc = 0.0004395 

^ 99 % 

SE q mc = 0.0013269 

v 99.9% 


5.4 cr = 25% 

In this section the value of the volatility is set to cr = 25%, while the other parameters are chosen as in 
Table EU 
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First let us display in Figure 20 the absolute and relative error of the variance, while in 14 are collected 
both the absolute and absolute value of relative error for mean and variance concerning the Vasicek interest 
rate model. 




Figure 20: Semilogy scale plot of the absolute error of the variance of R^ evaluated at time T = 1 (left) 
and the absolute relative error of the variance (right) of R^' 1 at time T = 1. In both cases the parameters 
of the Vasicek interest rate model are a = 0.1, /3 = 0.2, a = 0.25. 


Degree of PCE Average Error Variance Error Average relative error Variance relative error 


0 

8.2289e-08 

5.1512e-02 

9.1279e-10 

1 .0000e+00 

1 

8.2346e-08 

3.4693e-03 

9.1342e-10 

6.7349e-02 

2 

8.2320e-08 

3.4693e-03 

9.1314e-10 

6.7349e-02 

3 

8.2375e-08 

3.4693e-03 

9.1374e-10 

6.7349e-02 

4 

8.2274e-08 

3.4693e-03 

9.1262e-10 

6.7349e-02 

5 

8.2351e-08 

3.4693e-03 

9.1348e-10 

6.7349e-02 

6 

8.2342e-08 

3.4693e-03 

9.1338e-10 

6.7349e-02 

7 

8.2285e-08 

3.4693e-03 

9.1275e-10 

6.7349e-02 

8 

8.2322e-08 

3.4693e-03 

9.1316e-10 

6.7349e-02 

9 

8.2348e-08 

3.4693e-03 

9.1345e-10 

6.7349e-02 

10 

8.2306e-08 

3.4693e-03 

9.1298e-10 

6.7349e-02 

11 

8.2321e-08 

3.4693e-03 

9.1314e-10 

6.7349e-02 

12 

8.2323e-08 

3.4693e-03 

9.1317e-10 

6.7349e-02 

13 

8.2356e-08 

3.4693e-03 

9.1353e-10 

6.7349e-02 

14 

8.2346e-08 

3.4693e-03 

9.1343e-10 

6.7349e-02 

15 

8.2359e-08 

3.4693e-03 

9.1357e-10 

6.7349e-02 


Table 14: Absolute error of the average and the variance of at time T = 1, with respect to Vasicek 
parameters a = 0.1, f3 = 0.2, a = 25%. 


As for a = 15% the Table 14 shows the stationary behavior of the absolute and relative errors of mean 
and variance. Such features can be motivated as in the a = 25% case, indeed both the same computations 
and bounds hold again. 
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Coming back to the values in Table 11 let us compute a Monte Carlo sampling of size 5000 of 


( p) 


with p = 15, whose histogram, see fig. m is compared with the probability density function of the normal 
random variable Rt . Moreover in the same fig. [21] we show a sampling obtained exploiting the standard 
Monte Carlo technique. 




Figure 21: Left plot: Probability density function of the Vasicek interest rate model, with parameters 
a = 0.1, /3 = 0.2, cr = 0.25 and Rq = 110, (blue curve) and histogram of a Monte Carlo sampling (size 
= 5000) of the R.^ 1 for p = 15. With the same parameters, in the right plot we compare the analytical 
probability density function of the Vasicek model, with the one obtained using the Monte Carlo sampling 
of size 5000, for Rt , being T = 1. 


The errors of quantiles for 7 = 99%, resp. for 7 


99.9%, are shown in Fig. 22 and in Table 15 


The quantile computed by means of a standard Monte Carlo sampling of the analytical solution, whose 
size is equal to the one used for PCE, gives the following standard error 

SE 0 mc = 0.0008740 , SE q mc = 0.0022574 . 

^99% ^99.9% 
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Degree of PCE 

e 99% 

e 99.9% 

0 

1.1984e-01 

1.5919e-01 

1 

3.9090e-01 

5.2498e-01 

2 

3.9090e-01 

5.2498e-01 

3 

3.9090e-01 

5.2498e-01 

4 

3.9090e-01 

5.2498e-01 

5 

3.9090e-01 

5.2498e-01 

6 

3.9090e-01 

5.2498e-01 

7 

3.9090e-01 

5.2498e-01 

8 

3.9090e-01 

5.2498e-01 

9 

3.9090e-01 

5.2498e-01 

10 

3.9090e-01 

5.2498e-01 

11 

3.9090e-01 

5.2498e-01 

12 

3.9090e-01 

5.2498e-01 

13 

3.9090e-01 

5.2498e-01 

14 

3.9090e-01 

5.2498e-01 

15 

3.9090e-01 

5.2498e-01 


Table 15: Absolute errors of the two quantiles Qgg% and Qgg.g% of the PCE approximation of the Vasicek 
Interest rate model at time T = 1. The parameters are set as a = 0.1, /3 = 0.2, a = 0.25 and Rg = 110. 


Quantile absolute error 



Figure 22: Absolute errors of the two quantiles Qgg% and Qgg.g% of the PCE approximation of the Vasicek 
interest rate model at time T = 1. The parameters are set as a = 0.1, (3 = 0.2, a = 0.25 and Rg = 110. 


5.5 <r = 30% 

In this section the value of the volatility is set to a = 30%, while the other parameters are chosen as in 
Table E] 

First let us display in Figure [23] the absolute and relative error of the variance, while in Table [l6| are 
collected both the absolute and absolute value of relative error for mean and variance concerning the Vasicek 
interest rate model. 
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Figure 23: Semilogy scale plot of the absolute error of the variance of evaluated at time T = 1 (left) 
and the absolute relative error of the variance (right) of at time T = 1. In both cases the parameters 
of the Vasicek interest rate model are a = 0.1, f3 = 0.2, a = 30%. 


Degree of PCE Average Error Variance Error Average relative error Variance relative error 


0 

8.2289e-08 

7.4178e-02 

9.1279e-10 

1 .0000e+00 


1 

8.2400e-08 

4.9958e-03 

9.1402e-10 

6.7349e-02 


2 

8.2292e-08 

4.9958e-03 

9.1282e-10 

6.7349e-02 


3 

8.2276e-08 

4.9958e-03 

9.1264e-10 

6.7349e-02 


4 

8.2334e-08 

4.9958e-03 

9.1328e-10 

6.7349e-02 


5 

8.2374e-08 

4.9958e-03 

9.1374e-10 

6.7349e-02 


6 

8.2287e-08 

4.9958e-03 

9.1277e-10 

6.7349e-02 


7 

8.2332e-08 

4.9958e-03 

9.1326e-10 

6.7349e-02 


8 

8.2349e-08 

4.9958e-03 

9.1345e-10 

6.7349e-02 


9 

8.2326e-08 

4.9958e-03 

9.1320e-10 

6.7349e-02 


10 

8.2346e-08 

4.9958e-03 

9.1342e-10 

6.7349e-02 


11 

8.2296e-08 

4.9958e-03 

9.1286e-10 

6.7349e-02 


12 

8.2324e-08 

4.9958e-03 

9.1318e-10 

6.7349e-02 


13 

8.2349e-08 

4.9958e-03 

9.1346e-10 

6.7349e-02 


14 

8.2309e-08 

4.9958e-03 

9.1301e-10 

6.7349e-02 


15 

8.2345e-08 

4.9958e-03 

9.1342e-10 

6.7349e-02 


Table 16: Absolute error of the average and the variance of R^ at time T = 

1 , with respect to the 

Vasicek 

parameters a 

As for a = 

= 0.1, /3 = 0.2, cr = 30%. 

15% and a = 25%, the Table 16 shows the stationary behavior of the absolute and 

relative 

errors of mean and variance, which can 

be motivated for the same reasoning 

for a = 15%, and we 

recover 

the same computations and bounds obta 
Coming back to the values in Table 

iin( 

ii 

3 d for the previous values of a. 

(v) 

let us compute a Monte Carlo sampling of size 5000 of Rf , 


with p = 15, whose histogram, see fig. |24[ is compared with the probability density function of the normal 
random variable Rt■ Moreover in the same fig. [24] we show a standard Monte Carlo sampling. 
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Figure 24: Left plot: Probability density function of the Vasicek interest rate model, with parameters 
a = 0.1, /? = 0.2, cr = 0.30 and Rq = 110, (blue curve) and histogram of a Monte Carlo sampling (size 
= 5000) of the for p = 15. With the same parameters, in the right plot, we compare the analytical 
probability density function (24), with the one obtained using the Monte Carlo sampling of size 5000, for 
Rt , being T = 1 . 


The errors of quantiles for 7 = 99%, resp. for 7 


99.9%, are shown in Fig. 25 and in Table 17 


Degree of PCE 

e 99% 

£99.9% 

0 

1.7256e-01 

2.2923e-01 

1 

4.4032e-01 

5.9178e-01 

2 

4.4032e-01 

5.9178e-01 

3 

4.4032e-01 

5.9178e-01 

4 

4.4032e-01 

5.9178e-01 

5 

4.4032e-01 

5.9178e-01 

6 

4.4032e-01 

5.9178e-01 

7 

4.4032e-01 

5.9178e-01 

8 

4.4032e-01 

5.9178e-01 

9 

4.4032e-01 

5.9178e-01 

10 

4.4032e-01 

5.9178e-01 

11 

4.4032e-01 

5.9178e-01 

12 

4.4032e-01 

5.9178e-01 

13 

4.4032e-01 

5.9178e-01 

14 

4.4032e-01 

5.9178e-01 

15 

4.4032e-01 

5.9178e-01 


Table 17: Absolute errors of the two quantiles Qg$% and Qgg.g% of the PCE approximation of the Vasicek 
Interest rate model at time T = 1. The parameters are set as a = 0.1, /3 = 0.2, a = 0.30 and Rq = 110. 


The quantiles computed by means of a standard Monte Carlo sampling of the analytical solution, whose 
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Quantile absolute error 



Figure 25: Absolute errors of the two quantiles Q 99 % and Qgg. 9 % of the PCE approximation of the Vasicek 
interest rate model at time T = 1. The parameters are set as a = 0.1, j3 = 0.2, a = 0.30 and Rq = 110. 


size is equal to the one used for PCE, give the following standard error 

SE q mc = 0.0010324 

^99% 

SE 0 mc = 0.0025387 

^99.9% 


6 CIR interest rate model 


In what follows we consider the PCE method for the CIR interest rate model, introduced in 1985 by John 
C. Cox, Jonathan E. Ingersoll and Stephen A. Ross, see m , and defined by the following SDE 


dRt = (a — /3Rt) dt + ay/RtdWt 


t j 


(34) 


where a, /3 and er are positive, real-valued constants, and Wt is a standard Wiener process. Let us recall that 
eq. (34) generalizes the Vasicek model seen in Sec. [ 5 ] in fact, even if the drift term (a — f3Rt) is common 
between the twos, the standard deviation factor is different and, in particular, it equals a\fR~ t for the CIR 
model ensuring to avoid negative interest rates if a and /3 are positive. Moreover if 2a/3 > <r 2 , we cannot 
have Rt to be equal to zero. It is worth to mention that, unlike in the case of the Vasicek model, there is 
not an analytical solution to eq. ( [34| . The CIR process is ergodic and possesses a stationary distribution, 
see, e.g., m Section 3.2.3] and fToTSection 23.5] for further details. 

The statistics of the CIR model for fixed T > 0, read as follows 


E[R T ]=R(0)e-^ T +^(l-e-^ T ), 

Var [R t ] = jR( 0)(e“^ - e" 2 ^) + ^ (l - 2e^ T + e" 2 ^) . 


6.1 PCE approximation 

In order to apply the PCE method to the process Rt in ( |34| evaluated at a specific time T > 0, let us 
consider the stochastic process {V*}t>o 

Y t = y/Rt, (35) 
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then, applying the Ito-Doebling Lemma, Y t satisfies the following SDE 


dY t = fa - pY t 2 - -A dt + l -adW t , 


( 36 ) 


therefore we have canceled the interactions between the state of the original process R t and the increments 

Chapter 2] and [T5] Chapter 4] for 


of the Wiener process. The transformation defined by (351, see |14l 
details, allows to reduce simulations instability as well as it simplifies the application of Doss approach, 
developed in Sec. |3.2| Indeed the volatility is constant thus its derivative is null. 

Actually there is a subtleties concerning the drift p(x) = ^ [a — fix 2 — -ler 2 ), it is Lipschitz on intervals 
(a, +oo), where a > 0. Since the solution of eq. ( [36] ) is strictly positive, providing 2a/3 > a 2 , it is quite 
reasonable to apply the Doss theory in such setting. 


As seen studying the gBm model, see Sec. [4] and the Vasicek model, see Sec. [5j the solution to (351 can 
be expressed in terms of a functional H 


Y t = H(D t , W t ) 


(37) 


for every fixed time T, which is often referred as maturity time. In particular the scalar function H(x,y) is 
determined by solving the ODE (141, namely 

H(x,y) = ^V + x . 

Then the process {D t } t > 0 is path-wise determined by solving the following ODE 

1 


D{U}) 2 H(D{u),W t (u)) 

D(uj) = R 0 


(a~f3{H(D(cu),W t (oo))) 2 ~\a 2 ) t> 0 

t = 0 


(38) 


for almost all oj € ft. Notice that the dependence on t of the function D{uj), for each fixed sample event w, 
is omitted in order to shorten the notation. Moreover D(w) represents its derivative with respect to time. 
As in Sec. [4] and Sec. [5] the Wiener process is 


W t ~Af(0,t) ,t € [0,T], 


(39) 


which allows to set W t = y/2t£, and we still consider the basic random variable s £ ~ A/”(0,1/2), as made in 
the previous gBm and Vasicek sections, to fulfill the restrictions imposed by the software we have used to 
numerically implement the PCE method. In particular (381 is numerically solved, for a fixed sample event 
ui, by means of an adaptive Runge-Kutta method of order 4 (RK4), where the absolute tolerance, resp. 
relative tolerance, is set as le-7, resp. le-5, see, e.g., |361 Section II.1, Scetion II.4] for further details. 


We can therefore define the solution of (381 at time T as Dt = A4i(£, 0), where 0 £ R 3 collects all the 


parameters that characterized the dynamics of the CIR interest rate model (341. Then using both (37) and 

2 


([35l, we get 


0) :=R t = (^Wt + Dt) = (|V2T£ + AM£,e)) . 


(40) 


We shorten the notation omitting the explicit dependence of the CIR model on the parameters vector 
0 = (a,/3,cr) € K 3 , being a,/3,cr positive constants. Applying the truncated, at degree p , PCE to Rt, we 
obtain 

p 

?(p) 


=Y,Ci*i{0 , 


»=o 


where the coefficients are detecting using (10), evaluating Rt = A4(£) on a set of realizations of the basic 
random variable, namely we compute (401 at taking N = p and referring to these values as 


{Rtj }?=1 = {M&)}f= 1 

where such realizations are determined in two steps: 


(41) 
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each Gaussian quadrature nodes, belonging to can be interpreted as the image of a suitable 

ujj £ VL through £. Therefore by integrating (381 for each path in {wi,...,wat} C Q we achieve a 


set of realization of Dt, let us say {7W i ) }y=i - It is worth to mentioned that we are integrating N 

independent ODEs. 


By means of (40) we get the required set {Rtj }{Li, where the second summand is simply {^y/2T£j}^ =1 . 


6.2 Numerical application: overview of the computations 

Even if the CIR model does not admit close solution, see, e.g. mi Section 4.4.4], it possesses a stationary 
distribution, namely, by following m Section 3.2.3], Rt features a non-central chi-squared distribution. In 
particular let us denote by f^ T the analytical probability density function of Rt, for a positive end time T , 
then 

fR t(P^) Ix 2 (q,^t)/ct ^ ^ [0, +oo) , (42) 

where the subscript denotes the random variable considered. Moreover the constants are defined as 

= 4/3 

° T <t 2 (1 — exp(—/3T)) 

4a 

q = ~2 > 

\ T = c t Rq exp(—/3T) , 

where the last two constants are, respectively, the degree of freedom and the non-centrality parameter, see, 
e.g., m Section 3.2.3], for further details. In order to achieve a natural number q = 3 let us set a = |cr 2 , 
where a is considered as a known data. 


(43) 

(44) 

(45) 


Therefore to compute the probability density function and the analytic quantiles we set a in order to 


satisfy (441, where <7 = 3, moreover the PCE-approximation is computed for a set of volatility, namely 


a = {15%, 25%, 30%} 

while the other parameters are shown in Table |T8| 


Parameters 

P 

Rq 

T 

Values 

0.2 

110 

2 


Table 18: Parameters involved in the CIR interest model. 


Then for each value of a and for an increasing set of degree p £ {0,1,..., 15} the absolute errors of the 
mean, resp. of the variance, are computed as follows 


,(p) 

t MEAN 

Ip) 

t VAR 


IE [R7 1 ] — IE 


R 


(p) 


Var [Rt] — Var R 


jtp) 


furthermore the absolute value of the relative errors are considered, where, as we have assumed discussing 
the gBm and the Vasicek model, due to the logarithmic scale, we report their absolute values, to highlight 
its order, rather that its numerical value. These errors are given by 


RE 


(p) 

MEAN 


RE 


(p) 

VAR 


Jp) 

t MEAN 

E[I? T ] 

,(p) 

t V AR 

Var [I?t] 


Lastly we compute the two quantiles Q 7 , where 7 = 99% and 7 = 99.9%, of the PCE approximation 
R^\ where p = 0,1,..., 15. Proceeding as in Sec. [ 4 J we consider the 7 -th sample quantile, namely the 
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(K + l)-th realization of the sampling of R^\ sorted in ascending order, such that K < [ 7 M], where M is 
the size of the sampling, and [•] denotes the integer part of the real number within the brackets. In what 
follows we always employ a Latin Hypercube Sampling (LHS) technique of R?\ see |241 [25] and also [28] 
for further references, of size M = 5000. In particular we have A' g g% = 4951, while K 999 % = 4996, and we 


compute the absolute error of Q 1 versus the analytical values Q 1 . Therefore we evaluate e 7 = 


Q~f Q *i 


where Q 7 is the quantile of a normal random variable of mean E [Ay] and variance Var [Ay] and, indicating 
its cumulative density function by F} j T *, we have Q 7 = Fff^ ( 7 ). 


We compare these quantiles with the ones estimated using a standard Monte Carlo sampling, whose 
size is M = 5000, of the analytical solution Rt . In particular we compute the standard error for L = 200 
independent estimates {Qif c (/)} of the quantile. Thus 


se q 


MC 

99% 


Vl 


where a estimates the standard deviation of L independent estimates of Q 7 , thus 


1 

L- 1 


E K c (') - 


while Q h ' C 


represents the arithmetic mean of {Q(f c {l)} 


L 

l=i‘ 


Remark 7 . As for the gBm and the Vasicek case, the choice of the basic estimator for quantiles, i.e. the 
sample quantile, is due at focusing our attention to the efficacy of the method, instead that on the accuracy 
of the estimates as well as providing a fair comparison between the data achieved. More accurate techniques 
are, e.g., those presented in the L-estimator, and also the Harrel Davis (HD) estimators, see, e.g., 
for further details. 


6.3 a - 15% 

In this section the value of the volatility is set to a 
for q = 3. Its numerical value is 


15%, therefore as discussed in Sec 6.2 a satisfies (441, 


a = 0.005625 , 

while the other parameters are chosen as in Table [78] 


First let us display in Figure [26] Figure [27] and Table [T9] the absolute and relative error of the average, 
resp. variance, for the CIR interest rate model. 
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Figure 26: Semilogy scale plot of the absolute error of the mean (left) and the variance (right) computed 
via PCE for the CIR at time T = 2, whose parameters are a = 0.005625, /? = 0.2, a = 15% and starting 
value R 0 = 110, for a set of degrees p = {0,1,2,, 15}. 



Figure 27: Semilogy scale plot of the absolute value of the relative error of the mean (left) and the variance 
(right) computed via PCE for the CIR at time T = 2, whose parameters are a = 0.005625, /3 = 0.2, a — 15% 
and starting value Rq = 110, for a set of degrees p = (0,1,2,..., 15}. 


The absolute and relative errors of mean and variance are stationary. Therefore there is an error that 
corrupts the spectral convergence of PCE-approximation. In particular the numerical values suggest that 
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Degree of PCE Average Error Variance Error Average relative error Variance relative error 


0 

9.2721e-03 

2.7353e+00 

1.2572e-04 

1 .0000e+00 

1 

6.2470e-04 

1.8449e-01 

8.4691e-06 

6.7449e-02 

2 

6.2470e-04 

1.8434e-01 

8.4691e-06 

6.7395e-02 

3 

6.2470e-04 

1.8434e-01 

8.4691e-06 

6.7395e-02 

4 

6.2470e-04 

1.8434e-01 

8.4691e-06 

6.7395e-02 

5 

6.2470e-04 

1.8434e-01 

8.4691e-06 

6.7395e-02 

6 

6.2470e-04 

1.8434e-01 

8.4691e-06 

6.7395e-02 

7 

6.2470e-04 

1.8434e-01 

8.4691e-06 

6.7395e-02 

8 

6.2470e-04 

1.8434e-01 

8.4691e-06 

6.7395e-02 

9 

6.2470e-04 

1.8434e-01 

8.4691e-06 

6.7395e-02 

10 

6.2470e-04 

1.8434e-01 

8.4691e-06 

6.7395e-02 

11 

6.2470e-04 

1.8434e-01 

8.4691e-06 

6.7395e-02 

12 

6.2470e-04 

1.8434e-01 

8.4691e-06 

6.7395e-02 

13 

6.2470e-04 

1.8434e-01 

8.4691e-06 

6.7395e-02 

14 

6.2470e-04 

1.8434e-01 

8.4691e-06 

6.7395e-02 

15 

6.2470e-04 

1.8434e-01 

8.4691e-06 

6.7395e-02 


Table 19: Absolute error of the average and the variance of Rj^ at time T = 2, with respect to the CIR 
parameters a = 0.005625, /3 = 0.2, a = 15%. 


it is not due to numerical error made in order to compute {Rr,j}jL i> but it is due to the approximation of 
the solution (34) obtained by means of (|40| . 

Coming back to the values in Table |l8jlet us compute a sampling of size 5000 of R^ \ with p = 15, whose 
histogram, see fig. |28| is compared with the probability density function of the normal random variable Rt- 
Moreover in the same fig. [28]we show a sampling obtained exploiting the standard Monte Carlo technique. 




Figure 28: Left plot: Probability density function of the CIR interest rate model, with parameters a = 
0.005625, /3 = 0.2, a = 0.15 and Rq = 110, (blue curve) and histogram of a Monte Carlo sampling (size 
= 5000) of the R^ for p = 15. With the same parameters, in the right plot we compare the analytical 
probability density function of the CIR model, with the one obtained using the Monte Carlo sampling of 
size 5000, for Rt , being T = 2. 
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The errors of quantiles for 7 = 99%, resp. for 7 


99.9%, are shown in Fig. 29 and in Table 20 


Degree of PCE 

e 99% 

e 99.9% 

0 

3.9106e+00 

5.3240e+00 

1 

1.8048e-01 

3.3009e-01 

2 

1.4218e-01 

2.5449e-01 

3 

1.4218e-01 

2.5449e-01 

4 

1.4218e-01 

2.5449e-01 

5 

1.4218e-01 

2.5449e-01 

6 

1.4218e-01 

2.5449e-01 

7 

1.4218e-01 

2.5449e-01 

8 

1.4218e-01 

2.5449e-01 

9 

1.4218e-01 

2.5449e-01 

10 

1.4218e-01 

2.5449e-01 

11 

1.4218e-01 

2.5449e-01 

12 

1.4218e-01 

2.5449e-01 

13 

1.4218e-01 

2.5449e-01 

14 

1.4218e-01 

2.5449e-01 

15 

1.4218e-01 

2.5449e-01 


Table 20: Absolute errors of the two quantiles Qgg% and < 599 . 9 % of the PCE approximation of the CIR 
Interest rate model at time T = 2. The parameters are set as a = 0.005625, (3 = 0.2, cr = 0.15 and 
R 0 = 110. 


Quantile absolute error 



Figure 29: Absolute errors of the two quantiles <5gg% and <5gg.g% of the PCE approximation of the CIR 
interest rate model at time T = 2. The parameters are set as a = 0.005625, j3 = 0.2, a = 0.15 and Rq = 110. 


The quantiles computed by means of a standard Monte Carlo sampling of the analytical solution, whose 
size is equal to the one used for PCE, give the following standard error 


SE q mc = 0.0201547 

^99% 

SE q mc = 0.0534291 

^ 99 . 9 % 
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6.4 <r = 25% 


In this section the value of the volatility is set to a 
get 


25%, consequently by means of (441, where q = 3, we 


a = 0.046875 . 


The other parameters are chosen as in Table [18] First let us display in Figure [30] Figure [31] and Table [21] 
the absolute and relative error of the average, resp. variance, for the CIR interest rate model. 




Figure 30: Semilogy scale plot of the absolute error of the mean (left) and the variance (right) computed 
via PCE for the CIR at time T = 2, whose parameters are a = 0.046875, /I = 0.2, a = 25% and starting 
value i? 0 = HO, for a set of degrees p = {0,1, 2,..., 15}. 


Degree of PCE 

Average Error 

Variance Error 

Average relative error 

Variance relative error 

0 

2.5756e-02 

7.6005e+00 

3.4906e-04 

1 .0000e+00 

1 

1.7373e-03 

5.1400e-01 

2.3537e-05 

6.7627e-02 

2 

1.7373e-03 

5.1284e-01 

2.3537e-05 

6.7475e-02 

3 

1.7373e-03 

5.1284e-01 

2.3537e-05 

6.7475e-02 

4 

1.7373e-03 

5.1284e-01 

2.3537e-05 

6.7475e-02 

5 

1.7373e-03 

5.1284e-01 

2.3537e-05 

6.7475e-02 

6 

1.7373e-03 

5.1284e-01 

2.3537e-05 

6.7475e-02 

7 

1.7373e-03 

5.1284e-01 

2.3537e-05 

6.7475e-02 

8 

1.7373e-03 

5.1284e-01 

2.3537e-05 

6.7475e-02 

9 

1.7373e-03 

5.1284e-01 

2.3537e-05 

6.7475e-02 

10 

1.7373e-03 

5.1284e-01 

2.3537e-05 

6.7475e-02 

11 

1.7373e-03 

5.1284e-01 

2.3537e-05 

6.7475e-02 

12 

1.7373e-03 

5.1284e-01 

2.3537e-05 

6.7475e-02 

13 

1.7373e-03 

5.1284e-01 

2.3537e-05 

6.7475e-02 

14 

1.7373e-03 

5.1284e-01 

2.3537e-05 

6.7475e-02 

15 

1.7373e-03 

5.1284e-01 

2.3537e-05 

6.7475e-02 


Table 21: Absolute error of the average and the variance of R^' 1 at time T = 2, with respect to the CIR 
parameters a = 0.046875, /3 = 0.2, a = 25%. 


The absolute and relative errors of mean and variance are stationary, see Table |21| Therefore there is an 
error that corrupts the spectral convergence of PCE-approximation, and the obtained values suggest that it 
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Figure 31: Semilogy scale plot of the absolute value of the relative error of the mean (left) and the variance 
(right) computed via PCE for the CIR at time T = 2, whose parameters are a = 0.046875, /3 = 0.2, cr = 25% 
and starting value Rq = 110, for a set of degrees p = {0,1,2,..., 15}. 


is not due to numerical errors related to the computation of {Rr,j}j= i> but, instead, it is due to the choice 
of approximating the solution (34) by means of (401. 


Coming back to the values in Table 21l let us compute a sampling of size 5000 of R^ \ with p = 15, whose 


histogram, see fig. |32| is compared with the probability density function of the normal random variable Rt- 
Moreover in the same fig. 32 we show a sampling obtained exploiting the standard Monte Carlo technique. 


Monte Carlo sampling of PCE 



Standaid Monte-Carlo sampling of PCE 



Figure 32: Left plot: Probability density function of the CIR interest rate model, with parameters a = 0.046875, 
f3 = 0.2, a = 0.25 and Ra = 110, (blue curve) and histogram of a Monte Carlo sampling (size = 5000) of the R^ 
for p = 15. With the same parameters, in the right plot, we compare the analytical probability density function of 
the CIR model and the standard Monte Carlo sampling of size 5000, for Rt, being T = 2. 


44 

























































































































































The errors of quantiles for 7 = 99%, resp. for 7 


99.9%, are shown in Fig. 33 and in Table 22 


Degree of PCE 

e 99% 

e 99.9% 

0 

6.5644e+00 

8.8778e+00 

1 

3.3747e-01 

5.4443e-01 

2 

2.3108e-01 

3.3443e-01 

3 

2.3108e-01 

3.3443e-01 

4 

2.3108e-01 

3.3443e-01 

5 

2.3108e-01 

3.3443e-01 

6 

2.3108e-01 

3.3443e-01 

7 

2.3108e-01 

3.3443e-01 

8 

2.3108e-01 

3.3443e-01 

9 

2.3108e-01 

3.3443e-01 

10 

2.3108e-01 

3.3443e-01 

11 

2.3108e-01 

3.3443e-01 

12 

2.3108e-01 

3.3443e-01 

13 

2.3108e-01 

3.3443e-01 

14 

2.3108e-01 

3.3443e-01 

15 

2.3108e-01 

3.3443e-01 


Table 22: Absolute errors of the two quantiles Qgg% and Qgg.g% of the PCE approximation of the CIR 
Interest rate model at time T = 2. The parameters are set as a = 0.046875, (3 = 0.2, cr = 0.25 and 
R 0 = 110. 


Quantile absolute error 



Figure 33: Absolute errors of the two quantiles Qgg% and Qgg.g% of the PCE approximation of the CIR 
interest rate model at time T = 2. The parameters are set as a = 0.046875, j3 = 0.2, a = 0.25 and Rq = 110. 

The quantiles computed by means of a standard Monte Carlo sampling of the analytical solution, whose 
size equals the one used for PCE, gives SE n Mc = 0.0353593 and SE n Mc = 0.0913327. 

^ 99 % ^ 99 . 9 % 
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6.5 <r = 30% 


In this section the value of the volatility is set to a 
get 


30%, consequently by means of (441, where q 


a = 0.0675 , 


3, we 


while the other parameters are chosen as in Table [l8| First let us display in Figure [34} Figure [35] and Table 


23 the absolute and relative error of the average, resp. variance, for the CIR interest rate model. 




Figure 34: Semilogy scale plot of the absolute error of the mean (left) and the variance (right) computed 
via PCE for the CIR at time T = 2, whose parameters are a = 0.0675, /3 = 0.2, a = 30% and starting value 
R 0 = HO, for a set of degrees p = {0,1, 2,..., 15}. 


Degree of PCE Average Error Variance Error Average relative error Variance relative error 


0 

3.7089e-02 

1.0947e+01 

5.0250e-04 

1 .0000e+00 

1 

2.5035e-03 

7.4166e-01 

3.3903e-05 

6.7748e-02 

2 

2.5035e-03 

7.3927e-01 

3.3903e-05 

6.7530e-02 

3 

2.5035e-03 

7.3927e-01 

3.3903e-05 

6.7530e-02 

4 

2.5035e-03 

7.3927e-01 

3.3903e-05 

6.7530e-02 

5 

2.5035e-03 

7.3927e-01 

3.3903e-05 

6.7530e-02 

6 

2.5035e-03 

7.3927e-01 

3.3903e-05 

6.7530e-02 

7 

2.5035e-03 

7.3927e-01 

3.3903e-05 

6.7530e-02 

8 

2.5035e-03 

7.3927e-01 

3.3903e-05 

6.7530e-02 

9 

2.5035e-03 

7.3927e-01 

3.3903e-05 

6.7530e-02 

10 

2.5035e-03 

7.3927e-01 

3.3903e-05 

6.7530e-02 

11 

2.5035e-03 

7.3927e-01 

3.3903e-05 

6.7530e-02 

12 

2.5035e-03 

7.3927e-01 

3.3903e-05 

6.7530e-02 

13 

2.5035e-03 

7.3927e-01 

3.3903e-05 

6.7530e-02 

14 

2.5035e-03 

7.3927e-01 

3.3903e-05 

6.7530e-02 

15 

2.5035e-03 

7.3927e-01 

3.3903e-05 

6.7530e-02 


Table 23: Absolute error of the average and the variance of R^' 1 at time T = 2, with respect to the CIR 
parameters a = 0.0675, f) = 0.2, a = 30%. 

The absolute and relative errors of mean and variance are stationary. Therefore there is an error that 
corrupts the spectral convergence of PCE-approximation. In particular the values suggest that it is not due 
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Figure 35: Semilogy scale plot of the absolute value of the relative error of the mean (left) and the variance 
(right) computed via PCE for the CIR at time T = 2, whose parameters are a = 0.0675, (3 = 0.2, cr = 30% 
and starting value Rq = 110, for a set of degrees p = {0,1,2,..., 15}. 


to numerical error coming from the computations of { Rr.j}j _\, but, instead, it is due to the approximation 
of the solution (341 by means of (40). 


Coming back to the values in Table 18 let us compute a Monte Carlo sampling of size 5000 of R, 


(p) 
T > 


with p = 15, whose histogram, see fig. |36[ is compared with the probability density function of the normal 
random variable Rt■ Moreover in the same fig. [36] we show a standard Monte Carlo sampling. 


Monte Carlo sampling of PCE 



Standaid Monte-Carlo sampling of PCE 



Figure 36: Left plot: Probability density function of the CIR interest rate model, with parameters a = 0.0675, 
= 0.2, a = 0.30 and Ra = 110, (blue curve) and histogram of a Monte Carlo sampling (size = 5000) of the R^ 
for p = 15. With the same parameters, in the right plot we compare the analytical probability density function of 
the CIR model, with the one obtained using the Monte Carlo sampling of size 5000, for Rt, being T = 2. 
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The errors of quantiles for 7 = 99%, resp. for 7 


99.9%, are shown in Fig. 37 and in Table 24 


Degree of PCE 

e 99% 

e 99.9% 

0 

7.9085e+00 

1.0687e+01 

1 

4.3007e-01 

6.8095e-01 

2 

2.7688e-01 

3.7857e-01 

3 

2.7688e-01 

3.7857e-01 

4 

2.7688e-01 

3.7857e-01 

5 

2.7688e-01 

3.7857e-01 

6 

2.7688e-01 

3.7857e-01 

7 

2.7688e-01 

3.7857e-01 

8 

2.7688e-01 

3.7857e-01 

9 

2.7688e-01 

3.7857e-01 

10 

2.7688e-01 

3.7857e-01 

11 

2.7688e-01 

3.7857e-01 

12 

2.7688e-01 

3.7857e-01 

13 

2.7688e-01 

3.7857e-01 

14 

2.7688e-01 

3.7857e-01 

15 

2.7688e-01 

3.7857e-01 


Table 24: Absolute errors of the two quantiles Qgg% and < 599 . 9 % of the PCE approximation of the CIR 
Interest rate model at time T = 2. The parameters are set as a = 0.0675, (3 = 0.2, a = 0.30 and i? Q = 110. 


Quantile absolute error 



Figure 37: Absolute errors of the two quantiles <5gg% and <5g9.g% of the PCE approximation of the CIR 
interest rate model at time T = 2. The parameters are set as a = 0.0675, fl = 0.2, er = 0.30 and Rq = 110. 

The quantiles computed by means of a standard Monte Carlo sampling of the analytical solution, whose 
size is equal to the one used for PCE, give the following standard error 

SE n Mc = 0.0416182 

^99% 

SE 0 mc = 0.1128376 

V 99.9% 

Remark 8 . These computations, for each dynamics considered, point out that the volatility a influences 
the accuracy of PCE computations. Indeed the errors both of statistics and quantile increase proportionally 
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to a, since the randomness of the Brwonian motion becomes more relevant. 


Remark 9 In each case considered the lack of spectral convergence arises. Nevertheless such a point is 
not connected with the model nor with the related parameters, in fact it is a general consequence due to the 
approximations implied by the concrete application of the PCE-machinery used to compute the solution at 
a certain positive time T. Indeed the mean square convergence is ruled by the coefficients, see 0, therefore 
from a numerical point of view they are closed to machine precision, therefore after a certain degree they 
comes irrelevant. This motivates once more such behavior. 


7 Polynomial Chaos Expansion compared with MC and QMC 

In this section we compare the PCE method with the standard Monte Carlo (MC) method as well as 
with the guasTMonte Carlo (QMC) one. In particular we consider the accuracy in approximating the 
statistics of each model studied, namely the gBm equity model, the Vasicek and the CIR model, evaluated 
at time T, as well as the computational time costs to detect their mean and variance, as comparison criteria. 


It is worth to mention that both the MC and the QMC methods are usually applied to approximate 
analytical solutions, hence, from our interest models point of view we can consider them only for the gBm 
and the Vasicek case, while, for the CIR model, we regard them in relation with the Euler-Maruyama 
scheme. 

The usual Monte Carlo approximations require a set of M independent realizations of the 

process Xt- Then the mean and variance are determined as 


E [X T \ 


Var [X^] 


M 

me = jtY] x T ,i , 


l 


M ^' 

i =1 


_2 _ 
a MC ~ 


1 

M - 1 


M 

{Xr,i - Hmc ) 2 , 

i= 1 


(46) 

(47) 


then, exploiting the Law of large numbers , we have that both (461 and (471 , 
true value , hence their standard errors 


converges to the correspondent 




PMC 


m 


SE °hc = 


2 

MC 


n — 1 ’ 


are more informative than the single execution error. The QMC method uses low-discrepancy sequence 
{4>j}jLi generated by Sobol algorithm, that maximizes the uniformity of the sample points, see, e.g., [233 
Chapter 5], then it uses such values to approximate the statistics we are interested in, or, generally speaking, 
to compute integrals of the following type 


1 

g{x)f(x)dx « — g{4>i) , (48) 

1 i =0 

where f(x)dx is the measure induced by the considered random variable which characterizes the related 
image space. By setting g{x) = x, resp. g(x) = (x — p) 2 , we can approximate the mean of the random 
variable Xt, resp. its variance. Moreover, in QMC computations, the parameters for accuracy are given by 
the absolute errors 



^UQMC = |E [Xt] - pqmc I , 

Hmc = | Var [Xt] ~ Pqmc\ ■ 

We would like to underline that for both MC and QMC, we consider only the effective time to compute the 
statistics, without taking into account the amount of time spent for the detection of the grid of samplings 
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which is necessary to simulate the process. Moreover we highlight that in each of the case we consider, 
namely the gBm, resp. the Vasicek, resp. the CIR model, all the computations concerning the PCE method 
have been performed following the schemes proposed in sections [4] resp. [5] resp. [6j with related accuracy 
estimated by means of the absolute error of the average as well as of the variance. 


7.1 Geometric Brownian Motion 


Let us consider the SDE in eq. (16), where the parameters are taken as follows T = l,r = 3% while 


a = {15%, 25%, 30%} and So = 100. The MC and QMC methods approximate the statistics of the solution 
to the gBm model exploiting its analytical solution, see eq. (221, for an increasing size M of samplings 
M = {2 8 , 2 9 ,..., 2 16 } , and we then compute, for each value of M, the standard errors for the MC, resp. 
the absolute error for the QMC, are computed. In order to compare the accuracy and computational time 
costs of the three aforementioned methods, let us consider the plots in Fig. [38] Fig. [39] and Fig. [40] 
correspoinding to a = {15%, 25%, 30%}. 

Where the x-axis represents the computational time costs of the methods, while the error of the statis¬ 
tics, namely SE^ MC , SE^, e^ QMC , e CT ^ MC , ^mean and e VARi are displayed along the y-axis. Moreover 
the x-axis is normalized with respect to the highest data among MC, QMC and PCE values, actually on 
the x-axis we report related relative computational time costs is shown, which are more effective than the 
absolute ones. We also note that each point of the plot represents the computations for different number of 
realization points, namely M for MC and QMC, p for the PCE method. 




Figure 38: Semilogy scale plot which compares the accuracy and computational time costs for detecting the 
average (left) and variance (right) of the gBm by means of PCE, MC and QMC approaches. The parameters 
of the gBm equity model are r = 3%, a = 15%, T = 1 and Sq = 100. 
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Figure 39: Semilogy scale plot which compares the accuracy and computational time costs for detecting the 
average (left) and variance (right) of the gBm by means of PCE, MC and QMC approaches. The parameters 
of the gBm equity model are r = 3%, a = 25%, T = 1 and Sq = 100. 




Figure 40: Semilogy scale plot which compares the accuracy and computational time costs for detecting the 
average (left) and variance (right) of the gBm by means of PCE, MC and QMC approaches. The parameters 
of the gBm equity model are r = 3%, a = 30%, T — 1 and Sq = 100. 


Fig. [38] Fig. [39] and Fig. [40] point out the high accuracy as well as the very low computational effort 
required by PCE to get the solution. 

Let us display the values of computational time cost for the PCE, see Table [25] resp. for the MC and 
the QMC methods, see Table [26] as well as the error made using both the MC and the QMC methods, see 
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Table [27] 


of PCE p 

Time costs a = 15% 

Time costs cr = 25% 

Time costs cr = 

1 

0.0000 

0.0000 

0.0000 

2 

0.0000 

0.0160 

0.0000 

3 

0.0000 

0.0000 

0.0000 

4 

0.0160 

0.0000 

0.0160 

5 

0.0150 

0.0150 

0.0160 

6 

0.0150 

0.0160 

0.0320 

7 

0.0320 

0.0150 

0.0310 

8 

0.0310 

0.0150 

0.0310 

9 

0.0310 

0.0160 

0.0470 

10 

0.0150 

0.0310 

0.0320 

11 

0.0310 

0.0320 

0.0470 

12 

0.0310 

0.0310 

0.0310 

13 

0.0310 

0.0310 

0.0320 

14 

0.0460 

0.0310 

0.0310 

15 

0.0470 

0.0460 

0.0470 


Table 25: Time cost of computations for PCE, approximating gBm at T = 1 for r = 3%, 
a = {15%, 25%, 30%} and S 0 = 100. 


The computational time costs required for the detection of PCE’s coefficients, see eq. (101, and its 
statistics, see equations © and are irrelevant if compared with time required to evaluate the process, 
eq. (211, at {CjljLi, which are defined by Gaussian quadrature formula as stated in (101, this motivates 
the data in Table [25l 




a = 

15% 



a = 

25% 



a = 

30% 



Average 

Variance 

Average 

Variance 

Average 

Variance 

M 

MC 

QMC 

MC 

QMC 

MC 

QMC 

MC 

QMC 

MC 

QMC 

MC 

QMC 

256 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0.0150 

0.0000 

0.0000 

0.0160 

0.0160 

512 

0.0150 

0.0000 

0.0000 

0.0150 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0.0150 

1024 

0.0000 

0.0000 

0.0160 

0.0000 

0.0000 

0.0160 

0.0150 

0.0000 

0.0160 

0.0000 

0.0000 

0.0160 

2048 

0.0150 

0.0000 

0.0000 

0.0150 

0.0000 

0.0160 

0.0160 

0.0150 

0.0150 

0.0000 

0.0160 

0.0160 

4096 

0.0160 

0.0160 

0.0150 

0.0310 

0.0310 

0.0160 

0.0160 

0.0310 

0.0150 

0.0310 

0.0160 

0.0150 

8192 

0.0470 

0.0310 

0.0470 

0.0470 

0.0310 

0.0310 

0.0320 

0.0470 

0.0310 

0.0470 

0.0310 

0.0310 

16384 

0.0620 

0.0780 

0.0780 

0.0780 

0.0780 

0.0620 

0.0780 

0.0940 

0.0620 

0.0780 

0.0780 

0.0780 

32768 

0.1400 

0.1250 

0.1720 

0.1710 

0.1400 

0.1400 

0.1720 

0.1560 

0.1240 

0.1250 

0.1720 

0.1560 

65536 

0.2810 

0.2650 

0.3270 

0.3280 

0.2800 

0.2650 

0.3280 

0.3740 

0.2960 

0.2810 

0.3120 

0.3270 


Table 26: Computational time costs achieved by MC and QMC methods to detect the mean and variance 
of the process St, at T = 1, with parameters r = 3%, a = {15%, 25%, 30%} and So = 100. 
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Average 

Variance 


Size of sampling M 

se^ mc 

£ A»QMC 

SE a 2 

a MC 

e a 2 

a QMC 


256 

8.7186e-01 

1.7295e-01 

1.7234e+01 

6.5079e+00 


512 

6.9261e-01 

9.4193e-02 

1.5365e+01 

3.6819e+00 


1024 

4.9055e-01 

5.0744e-02 

1.0896e+01 

2.0760e+00 

cr = 15% 

2048 

3.4608e-01 

2.7110e-02 

7.6670e+00 

1.1656e+00 

4096 

2.4726e-01 

1.4388e-02 

5.5344e+00 

6.5155e-01 


8192 

1.6965e-01 

7.5954e-03 

3.6842e+00 

3.6254e-01 


16384 

1.2111e-01 

3.9921e-03 

2.6552e+00 

2.0082e-01 


32768 

8.5602e-02 

2.0905e-03 

1.8759e+00 

1.1076e-01 


65536 

6.1000e-02 

1.0913e-03 

1.3472e+00 

6.0845e-02 


256 

1.4618e+00 

3.1675e-01 

4.8443e+01 

2.5888e+01 


512 

1.1236e+00 

1.7307e-01 

4.0438e+01 

1.5299e+01 


1024 

8.1200e-01 

9.3620e-02 

2.9853e+01 

8.9624e+00 

cr = 25% 

2048 

5.7963e-01 

5.0254e-02 

2.1507e+01 

5.2056e+00 

4096 

4.1155e-01 

2.6811e-02 

1.5332e+01 

2.9990e+00 


8192 

2.8569e-01 

1.4233e-02 

1.0448e+01 

1.7149e+00 


16384 

2.0385e-01 

7.5252e-03 

7.5225e+00 

9.7390e-01 


32768 

1.4416e-01 

3.9648e-03 

5.3204e+00 

5.4964e-01 


65536 

1.0226e-01 

2.0828e-03 

3.7858e+00 

3.0844e-01 


256 

1.9748e+00 

3.9966e-01 

8.8419e+01 

4.4469e+01 


512 

1.4752e+00 

2.1905e-01 

6.9708e+01 

2.6743e+01 


1024 

9.5791e-01 

1.1890e-01 

4.1545e+01 

1.5910e+01 

cr = 30% 

2048 

6.8569e-01 

6.4066e-02 

3.0099e+01 

9.3693e+00 

4096 

4.9136e-01 

3.4316e-02 

2.1855e+01 

5.4664e+00 


8192 

3.4695e-01 

1.8292e-02 

1.5409e+01 

3.1625e+00 


16384 

2.4803e-01 

9.7119e-03 

1.1136e+01 

1.8158e+00 


32768 

1.7542e-01 

5.1389e-03 

7.8780e+00 

1.0354e+00 


65536 

1.2311e-01 

2.7113e-03 

5.4873e+00 

5.8685e-01 


Table 27: Errors of average and variance computed by MC and QMC methods for the gBm evaluated at 
time T = 1, with parameters r = 3%, a = {15%, 25%, 30%} and S 0 = 100. 


7.2 Vasicek interest rate model 


Proceeding as in Sec. 7.1 we consider the SDE (|23| characterizing the Vasicek model for the following set of 

10 -5 , a = {15%, 25%, 30%}, and we apply both the MC and QMC 
(|23} , and to compute the statistics of Rt at time 
We use the same type of plot exploited in the 


values for its parameters a = 0.1, /3 = 2 • 
method to approximate the analytical solution (241 to eq. 

T = 1 for an increasing size of samplings M = } 2 s ,..., 2 16 j 
previous section, in order to show the accuracy, the computational time costs and the errors of the PCE, 
resp. the MC , resp. the QMC, method. 

In particular Fig. [44] Fig. [42] and Fig. [43] point out the high accuracy of PCE and the corresponding 
low computational effort to get these results: the key points are the low number of simulation required to 


get the solution of eq. (231. 
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Figure 41: Semilogy scale plot which compares the accuracy and computational time costs for detecting 
the average (left) and variance (right) of the Vasicek interest rate model by means of PCE, MC and QMC 
approaches. The parameters are a = 0.1, /? = 2 ■ 10 —5 , a = 15%, T = 1 and i? 0 = 110. 




Figure 42: Semilogy scale plot which compares the accuracy and computational time costs for detecting 
the average (left) and variance (right) of the Vasicek interest rate model by means of PCE, MC and QMC 
approaches. The parameters are a = 0.1, (3 = 2- 10 -5 , a = 25%, T = 1 and Rq = 110. 


In order to underline the efficiency and accuracy of the PCE method, let us show the related compu¬ 
tational time costs, see Table [28] compared to those of the MC and QMC methods, see Table [30] also we 
highlight the numerical values concerning the absolute errors achieved using PCE method, see Table [29] as 
well as MC and QMC approaches, see Table [31] 
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Figure 43: Semilogy scale plot which compares the accuracy and computational time costs for detecting 
the average (left) and variance (right) of the Vasicek interest rate model by means of PCE, MC and QMC 
approaches. The parameters are a = 0.1, /? = 2 • 10 —5 , a = 30%, T = 1 and i? 0 = 110. 


As observed in the gBm-case the computational time costs required for the detection of PCE’s coefficients, 
equation (101, and its statistics, equ atio ns ([8]) and ([£]), are irrelevant if compared with time required to get 
the evaluation of the process in eq. (281 at quadrature nodes {£j}jLi- 


of PCE p 

Time costs a = 15% 

Time costs a = 25% 

Time costs cr 

0 

0.0030 

0.0030 

0.0030 

1 

0.0050 

0.0040 

0.0050 

2 

0.0080 

0.0060 

0.0080 

3 

0.0100 

0.0080 

0.0100 

4 

0.0140 

0.0100 

0.0140 

5 

0.0140 

0.0110 

0.0140 

6 

0.0160 

0.0150 

0.0160 

7 

0.0180 

0.0150 

0.0180 

8 

0.0200 

0.0160 

0.0200 

9 

0.0240 

0.0180 

0.0240 

10 

0.0270 

0.0190 

0.0270 

11 

0.0220 

0.0210 

0.0220 

12 

0.0220 

0.0240 

0.0220 

13 

0.0260 

0.0260 

0.0260 

14 

0.0260 

0.0270 

0.0260 

15 

0.0310 

0.0310 

0.0310 


Table 28: Time cost of computations for PCE, approximating the Vasicek interest rate model at T = 1 for 
a = 0.1,/? = 2 - 10- 5 ,<J = {15%, 25%, 30%} and R 0 = 110. 
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a = 

15% 

a = 

25% 

a = 

30% 

Degree of PCE p 

~jp) 

t MEAN 

~JF) 

k VAR 

tMEAN 

~J p) 

tVAR 

~Tp) 

tMEAN 

~Jp) 

tVAR 

0 

1.0006e-07 

2.2500e-02 

1.0006e-07 

6.2499e-02 

1.0006e-07 

8.9998e-02 

1 

1.0006e-07 

1.4999e-07 

1.0006e-07 

4.1664e-07 

1.0006e-07 

5.9996e-07 

2 

1.0006e-07 

1.4999e-07 

1.0006e-07 

4.1664e-07 

1.0006e-07 

5.9997e-07 

3 

1.0006e-07 

1.4999e-07 

1.0006e-07 

4.1664e-07 

1.0006e-07 

5.9996e-07 

4 

1.0006e-07 

1.4999e-07 

1.0006e-07 

4.1664e-07 

1.0006e-07 

5.9997e-07 

5 

1.0006e-07 

1.4999e-07 

1.0006e-07 

4.1664e-07 

1.0006e-07 

5.9996e-07 

6 

1.0006e-07 

1.4999e-07 

1.0006e-07 

4.1664e-07 

1.0006e-07 

5.9997e-07 

7 

1.0006e-07 

1.4999e-07 

1.0006e-07 

4.1664e-07 

1.0006e-07 

5.9996e-07 

8 

1.0006e-07 

1.4999e-07 

1.0006e-07 

4.1664e-07 

1.0006e-07 

5.9997e-07 

9 

1.0006e-07 

1.4999e-07 

1.0006e-07 

4.1664e-07 

1.0006e-07 

5.9996e-07 

10 

1.0006e-07 

1.4999e-07 

1.0006e-07 

4.1664e-07 

1.0006e-07 

5.9997e-07 

11 

1.0006e-07 

1.4999e-07 

1.0006e-07 

4.1664e-07 

1.0006e-07 

5.9996e-07 

12 

1.0006e-07 

1.4999e-07 

1.0006e-07 

4.1664e-07 

1.0006e-07 

5.9997e-07 

13 

1.0006e-07 

1.4999e-07 

1.0006e-07 

4.1664e-07 

1.0006e-07 

5.9997e-07 

14 

1.0006e-07 

1.4999e-07 

1.0006e-07 

4.1664e-07 

1.0006e-07 

5.9997e-07 

15 

1.0006e-07 

1.4999e-07 

1.0006e-07 

4.1664e-07 

1.0006e-07 

5.9997e-07 


Table 29: Absolut errors of PCE-approximation of the Vasicek interest rate model at T = 1 for a = 0.1, 
P = 2-10 ~ 5 ,a = {15%, 25%, 30%} and R 0 = 110. 




a = 

15% 



a = 

25% 



a = 

30% 



Average 

Variance 

Average 

Variance 

Average 

Variance 

M 

MC 

QMC 

MC 

QMC 

MC 

QMC 

MC 

QMC 

MC 

QMC 

MC 

QMC 

256 

0.0020 

0.0020 

0.0020 

0.0020 

0.0020 

0.0020 

0.0020 

0.0030 

0.0020 

0.0020 

0.0010 

0.0020 

512 

0.0020 

0.0030 

0.0010 

0.0020 

0.0020 

0.0020 

0.0010 

0.0030 

0.0020 

0.0030 

0.0030 

0.0020 

1024 

0.0030 

0.0030 

0.0020 

0.0030 

0.0040 

0.0050 

0.0030 

0.0030 

0.0030 

0.0030 

0.0030 

0.0040 

2048 

0.0050 

0.0060 

0.0050 

0.0060 

0.0070 

0.0050 

0.0050 

0.0060 

0.0040 

0.0050 

0.0050 

0.0050 

4096 

0.0080 

0.0110 

0.0100 

0.0130 

0.0110 

0.0100 

0.0110 

0.0100 

0.0080 

0.0080 

0.0100 

0.0090 

8192 

0.0160 

0.0170 

0.0210 

0.0210 

0.0170 

0.0170 

0.0200 

0.0230 

0.0180 

0.0170 

0.0220 

0.0220 

16384 

0.0290 

0.0350 

0.0350 

0.0400 

0.0340 

0.0340 

0.0380 

0.0430 

0.0330 

0.0310 

0.0400 

0.0390 

32768 

0.0580 

0.0640 

0.0740 

0.0820 

0.0640 

0.0700 

0.0790 

0.0890 

0.0630 

0.0660 

0.0770 

0.0820 

65536 

0.1220 

0.1170 

0.1570 

0.1430 

0.1260 

0.1300 

0.1560 

0.1580 

0.1170 

0.1380 

0.1410 

0.1560 


Table 30: Computational time costs achieved by MC and QMC methods to detect the mean and variance 
of the process Rt ,at T = 1 for a = 0.1, p = 2 • 10” 5 , a = {15%, 25%, 30%} and R 0 = 110. 
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Average 

Variance 


Size of sampling M 

se» mc 

£ HQMC 

SE a 2 

a MC 

€ (J 2 
° QMC 


256 

9.6455e-03 

1.4768e-03 

2.1093e-03 

3.6049e-04 


512 

6.9057e-03 

8.0728e-04 

1.5275e-03 

1.8175e-04 


1024 

4.6269e-03 

4.3579e-04 

9.6929e-04 

9.1525e-05 

cr = 15% 

2048 

3.3442e-03 

2.3302e-04 

7.1593e-04 

4.6045e-05 

4096 

2.3732e-03 

1.2367e-04 

5.0984e-04 

2.3146e-05 


8192 

1.6639e-03 

6.5245e-05 

3.5439e-04 

1.1627e-05 


16384 

1.1707e-03 

3.4252e-05 

2.4810e-04 

5.8373e-06 


32768 

8.2966e-04 

1.7907e-05 

1.7622e-04 

2.9291e-06 


65536 

5.8410e-04 

9.3297e-06 

1.2352e-04 

1.4692e-06 


256 

1.5392e-02 

2.4614e-03 

5.3712e-03 

1.0014e-03 


512 

1.1543e-02 

1.3455e-03 

4.2681e-03 

5.0487e-04 


1024 

7.9822e-03 

7.2632e-04 

2.8848e-03 

2.5424e-04 

cr = 25% 

2048 

5.5877e-03 

3.8837e-04 

1.9987e-03 

1.2790e-04 

4096 

3.9699e-03 

2.0612e-04 

1.4266e-03 

6.4294e-05 


8192 

2.7871e-03 

1.0874e-04 

9.9438e-04 

3.2297e-05 


16384 

1.9696e-03 

5.7086e-05 

7.0224e-04 

1.6215e-05 


32768 

1.3777e-03 

2.9845e-05 

4.8592e-04 

8.1365e-06 


65536 

9.7726e-04 

1.5549e-05 

3.4576e-04 

4.0811e-06 


256 

1.8935e-02 

2.9537e-03 

8.1284e-03 

1.4420e-03 


512 

1.3966e-02 

1.6146e-03 

6.2476e-03 

7.2701e-04 


1024 

9.3332e-03 

8.7158e-04 

3.9441e-03 

3.6610e-04 

cr = 30% 

2048 

6.5312e-03 

4.6604e-04 

2.7307e-03 

1.8418e-04 

4096 

4.6692e-03 

2.4735e-04 

1.9735e-03 

9.2584e-05 


8192 

3.3352e-03 

1.3049e-04 

1.4239e-03 

4.6508e-05 


16384 

2.3304e-03 

6.8504e-05 

9.8311e-04 

2.3349e-05 


32768 

1.6527e-03 

3.5815e-05 

6.9927e-04 

1.1717e-05 


65536 

1.1749e-03 

1.8659e-05 

4.9977e-04 

5.8768e-06 


Table 31: Errors of the average and variance computed by MC and QMC for the Vasicek interest rate model 
R t at time T = 1, with parameters a = 0.1, /3 = 2 • 10~' 5 , , a = {15%, 25%, 30%} and and R 0 = 110. 
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7.3 CIR interest rate model 


Even if the CIR model does not have a solution in closed form, we can still apply both the MC and the 
QMC methods, but considering the numerical solution of the SDE (34). 


As described in m, the Euler-Maruyama scheme applied to the transformed process SDE ( [36] ) is actually 
the Milstein scheme applied to the original SDE ( [34] ) . Hence let us integrate numerically the transformed 
SDE from 0 up to T = 2 using 200 time step. The mean and the variance are estimated, via MC and QMC, 
using a set of samplings of increasing size M = {2 8 ,..., 2 15 }, setting the CIR model parameters as follows 


/3 = 0.002, a = {15%, 25%, 30%}. Moreover exploiting results stated in Sec. 6.2 a = (3/4) • er 2 , thus we have 


a = {0.005625,0.046875,0.0675}. We note that, while the MC method applies straightforward, the QMC 
implementation requires some further discussion which is, in fact, close to the analysis made in section [XT] 
namely the numerical scheme allows to write the solution at T = 2 as a functional of all the independent 
increments of the Wiener process, one for each time step. 

Since the dimensionality of the problem equals the number of time steps, i.e. 200, then a multidimensional 
low-discrepancy sequence is computed, where the size M represents the number of points belonging to 
this sequence. Once the the SDE is solved the Euler-Maruyama scheme prescribes that each independent 
increment has to be replaced by an entry of a point belonging to such a multidimensional low-discrepancy 


sequence. Eventually, by applying (481 to the detected solution, we get the required statistics. 

in Fig. [44[ Fig. [45] and Fig. [46] we compare the computational time 


As we made in sections |7. 1 1 and [772 
costs and its related accuracy reached by the PCE method, resp. the MC and QMC methods. The data 
used to apply the PCE technique are computed following the approach developed in section [6[ We recall 
that the plot is characterized by setting along the x-axis the computational time costs of the methods, while 


the error of the statistics SE, 


^MC 5 


SE„ 


and - 


and the PCE absolute errors, are placed along 


1 MC ‘ u QMC 1 

the y-axis. Moreover the time costs are normalized with respect to the highest time value among MC, QMC 
and PCE, thus we display their relative computational time costs, hence gaining more informations than 
merely using the related absolute values. 




Figure 44: Semilogy scale plot which compare the accuracy and computational time costs for detecting the 
average (left) and variance (right) of the CIR model at T = 2 by means of PCE, MC and QMC approaches, 
with parameters a = 0.005625, /3 = 0.002 and a = 15%. 

Let us underline that Fig. [44] Fig. [45] and [46[ point out the high accuracy as well as the very low com¬ 
putational effort required by the PCE method to get the solution. We can appreciate such an effectiveness 
property of the PCE approach looking at the following tables, where the values of computational time cost 
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Figure 45: Semilogy scale plot which compare the accuracy and computational time costs for detecting the 
average (left) and variance (right) of the CIR model at T = 2 by means of PCE, MC and QMC approaches, 
with parameters a = 0.046875, /3 = 0.002 and a = 25%. 



Figure 46: Semilogy scale plot which compare the accuracy and computational time costs for detecting the 
average (left) and variance (right) of the CIR model at T = 2 by means of PCE, MC and QMC approaches, 
with parameters a = 0.0675, /3 = 0.002 and a = 30%. 


of PCE, see Table [32] MC and QMC, see Table [34| as well as the errors characterizing the PCE, see Table 
[33] MC and the QMC methods, see Table [35] 
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Degree of PCE p 

Time costs a = 15% 

Time costs a = 25% 

Time costs a = 30% 

0 

0.0030 

0.0050 

0.0040 

1 

0.0070 

0.0110 

0.0070 

2 

0.0100 

0.0090 

0.0100 

3 

0.0110 

0.0110 

0.0190 

4 

0.0150 

0.0210 

0.0160 

5 

0.0190 

0.0200 

0.0180 

6 

0.0190 

0.0210 

0.0190 

7 

0.0230 

0.0250 

0.0220 

8 

0.0260 

0.0290 

0.0260 

9 

0.0270 

0.0260 

0.0270 

10 

0.0320 

0.0290 

0.0310 

11 

0.0310 

0.0340 

0.0340 

12 

0.0390 

0.0380 

0.0360 

13 

0.0360 

0.0370 

0.0410 

14 

0.0410 

0.0410 

0.0430 

15 

0.0430 

0.0430 

0.0440 


Table 32: Computational time costs and absolute error of mean and variance for PCE-approximation of 
the CIR interest rate model evaluated at T = 2, whose parameters are a = {0.005625,0.046875,0.0675}, 
/3 = 0.002, cr = {15%, 25%, 30%} and R 0 = 110. 



a = 

15% 

a = 

25% 

a = 

30% 

Degree of PCE p 

~ip) 

t MEAN 

~Jp) 

t VAR 

~jp) 

t MEAN 

~Jp) 

t V AR 

~jp) 

t MEAN 

~Jp) 

t VAR 

0 

1.1227e-02 

4.9211e+00 

3.1187e-02 

1.3674e+01 

4.4910e-02 

1.9694e+01 

1 

7.6695e-06 

3.7012e-03 

2.3551e-05 

1.2358e-02 

3.5874e-05 

1.9850e-02 

2 

7.6696e-06 

3.4494e-03 

2.3552e-05 

1.0415e-02 

3.5876e-05 

1.5822e-02 

3 

7.6695e-06 

3.4494e-03 

2.3552e-05 

1.0415e-02 

3.5876e-05 

1.5822e-02 

4 

7.6696e-06 

3.4494e-03 

2.3552e-05 

1.0415e-02 

3.5876e-05 

1.5822e-02 

5 

7.6695e-06 

3.4494e-03 

2.3552e-05 

1.0415e-02 

3.5876e-05 

1.5822e-02 

6 

7.6696e-06 

3.4494e-03 

2.3552e-05 

1.0415e-02 

3.5876e-05 

1.5822e-02 

7 

7.6695e-06 

3.4494e-03 

2.3552e-05 

1.0415e-02 

3.5876e-05 

1.5822e-02 

8 

7.6696e-06 

3.4494e-03 

2.3552e-05 

1.0415e-02 

3.5876e-05 

1.5822e-02 

9 

7.6695e-06 

3.4494e-03 

2.3552e-05 

1.0415e-02 

3.5876e-05 

1.5822e-02 

10 

7.6696e-06 

3.4494e-03 

2.3552e-05 

1.0415e-02 

3.5876e-05 

1.5822e-02 

11 

7.6695e-06 

3.4494e-03 

2.3552e-05 

1.0415e-02 

3.5876e-05 

1.5822e-02 

12 

7.6696e-06 

3.4494e-03 

2.3552e-05 

1.0415e-02 

3.5876e-05 

1.5822e-02 

13 

7.6696e-06 

3.4494e-03 

2.3552e-05 

1.0415e-02 

3.5876e-05 

1.5822e-02 

14 

7.6696e-06 

3.4494e-03 

2.3552e-05 

1.0415e-02 

3.5876e-05 

1.5822e-02 

15 

7.6695e-06 

3.4494e-03 

2.3552e-05 

1.0415e-02 

3.5876e-05 

1.5822e-02 


Table 33: Absolute errors of PCE-approximation of the CIR interest rate model at T = 2 for 
a = {0.005625,0.046875,0.0675}, /3 = 0.002, a = {15%, 25%, 30%} and R 0 = 110. 


The computational time costs required for the detection of PCE’s coefficients, equation (101, and its 
statistics, equa tion s © and ©, are irrelevant if compared with time required to get the e valu ation of the 
process in eq. (40) at the quadrature nodes {£j}:?Li, this motivates the time data in Table 
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a = 

15% 



a = 

25% 



a = 

30% 



Average 

Variance 

Average 

Variance 

Average 

Variance 

M 

MC 

QMC 

MC 

QMC 

MC 

QMC 

MC 

QMC 

MC 

QMC 

MC 

QMC 

256 

0.0140 

0.0070 

0.0080 

0.0050 

0.0080 

0.0110 

0.0050 

0.0050 

0.0070 

0.0090 

0.0060 

0.0060 

512 

0.0090 

0.0090 

0.0070 

0.0060 

0.0080 

0.0110 

0.0070 

0.0070 

0.0090 

0.0080 

0.0080 

0.0060 

1024 

0.0110 

0.0120 

0.0110 

0.0100 

0.0120 

0.0150 

0.0100 

0.0110 

0.0120 

0.0140 

0.0100 

0.0100 

2048 

0.0190 

0.0200 

0.0190 

0.0170 

0.0230 

0.0240 

0.0180 

0.0210 

0.0200 

0.0200 

0.0180 

0.0170 

4096 

0.0330 

0.0390 

0.0340 

0.0320 

0.0370 

0.0370 

0.0320 

0.0310 

0.0380 

0.0380 

0.0360 

0.0360 

8192 

0.0750 

0.0670 

0.0620 

0.0590 

0.0740 

0.0790 

0.0650 

0.0580 

0.0790 

0.0680 

0.0590 

0.0610 

16384 

0.1410 

0.1500 

0.1190 

0.1980 

0.1330 

0.1370 

0.1260 

0.1170 

0.1570 

0.1660 

0.2010 

0.1450 

32768 

0.2610 

0.2690 

0.2350 

0.2390 

0.2630 

0.2810 

0.2540 

0.2580 

0.2840 

0.2690 

0.2410 

0.2430 


Table 34: Computational cost of computation of statistics of the CIR interest rate model Rt at T = 2, 
determined by MC, QMC, with a = {0.005625,0.046875,0.0675}, /3 = 0.002, a = {15%, 25%, 30%} and 
R 0 = 110. 




Average 

Variance 


Size of sampling M 

se» mo 

^MQMC 

SE a 2 

° MC 

e a 2 
° QMC 


256 

1.3515e-01 

1.5744e-02 

4.1411e-01 

3.5681e-01 


512 

9.8423e-02 

1.4373e-03 

3.1029e-01 

2.6892e-01 


1024 

6.9759e-02 

1.1029e-03 

2.2033e-01 

3.3319e-01 

a = 15% 

2048 

4.9068e-02 

1.3868e-04 

1.5413e-01 

1.6354e-01 

4096 

3.4559e-02 

3.7763e-04 

1.0811e-01 

7.1762e-02 


8192 

2.4551e-02 

4.2046e-04 

7.7155e-02 

8.7361e-02 


16384 

1.7269e-02 

1.8063e-04 

5.3984e-02 

8.6968e-03 


32768 

1.2218e-02 

1.2084e-04 

3.8214e-02 

4.4496e-02 


256 

2.3429e-01 

2.5344e-02 

1.2445e+00 

9.9417e-01 


512 

1.7451e-01 

1.7405e-03 

9.7546e-01 

7.6735e-01 


1024 

1.1406e-01 

1.0069e-03 

5.8900e-01 

9.3635e-01 

cr = 25% 

2048 

8.0223e-02 

6.4006e-04 

4.1199e-01 

4.5978e-01 

4096 

5.8942e-02 

8.0988e-04 

3.1449e-01 

2.0169e-01 


8192 

4.0456e-02 

4.7886e-04 

2.0951e-01 

2.4117e-01 


16384 

2.9153e-02 

3.2369e-04 

1.5385e-01 

2.2580e-02 


32768 

2.0581e-02 

3.1324e-04 

1.0844e-01 

1.2319e-01 


256 

2.6173e-01 

2.9876e-02 

1.5530e+00 

1.4338e+00 


512 

2.0446e-01 

1.6963e-03 

1.3391e+00 

1.1198e+00 


1024 

1.3572e-01 

7.0999e-04 

8.3396e-01 

1.3563e+00 

cr = 30% 

2048 

9.7037e-02 

1.0130e-03 

6.0278e-01 

6.6614e-01 

4096 

7.0453e-02 

1.0797e-03 

4.4931e-01 

2.9220e-01 


8192 

4.8704e-02 

4.4112e-04 

3.0364e-01 

3.4621e-01 


16384 

3.4759e-02 

4.0243e-04 

2.1871e-01 

3.1364e-02 


32768 

2.4443e-02 

4.4340e-04 

1.5295e-01 

1.7711e-01 


Table 35: Errors of the average and variance computed by MC and QMC for the CIR interest rate model 
R t at time T = 2, with parameters a = {0.005625,0.046875,0.0675}, /? = 0.002, ,a = {15%, 25%, 30%} 
and and R 0 = 110. 
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8 Conclusions 


In this paper we show how the Polynomial Chaos Expansion (PCE) method can be effectively used to 
approximate the solution of a given SDE, exploiting only some basic properties of the Brownian motion, 
which is assumed to drive the stochastic process we are interested in, and a suitable transformation of 
the original SDE we have taken into consideration. In particular we have applied the PCE technique to 
approximate the solutions of some of the most relevant interest rate models, namely the geometric Brownian 
motion model, the Vasicek model and the Cox-IngersolRRoss (CIR) model. 

The provided numerical examples deeply discuss the convergence properties of the PCE method in var¬ 
ious scenarios of volatility. Indeed detailed analysis on error properties of the approximated statistics is 
given when analytical solution is available also providing a comparison with respect to the related analytical 
probability density function. The latter is not the case for the CIR model, conversely to what happens for 
both the gBm and the Vasicek models. Hence, in the CIR case, the PCE approach is even more interesting 
also because it gives a general method to study those stochastic processes which are defined by SDEs which 
do not have solution in closed form. 

We show by concrete examples how the PCE method overcomes the performances of the standard Monte 
Carlo (MC) method as well as those of the quasi -Monte Carlo (QMC) one. In particular a close comparison 
between PCE, MC and QMC techniques, points out the advantage in using the PCE approach, especially 
in terms of computational costs. The latter result is due to the fact hat the PCE method requires just few 
realization of the random variable we want to approximate, therefore, even if the basic operations required 
by the PCE approach, are rather more time consuming, than those implied by the MC or the QMC ap¬ 
proach ( weighted average versus a quadrature formulas), the computational advantages in using the PCE 
technique are clear as shown, e.g. , by the CIR model analysis. 

Last but not least, it is worth to mention that the PCE approach can be split in two distinct section: 
the off-line computations and the real PCE approximation of the random variable we want to study. Thus 
the basic random variable, the orthogonal polynomials and some data used in the Gaussian quadrature 
formulas, are computed once and for all, since they do not depend on the particular model, e.g. the gBm, 
the Vasicek or the CIR modelm as we made in our analysis. Hence a clear saving concerning the required 
computational efforts is achieved every time the model is sensible to time variations, as well as in the case 
when a new calibration of the involved parameters is required. 
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